Decay of a superfluid currents in a moving system of strongly interacting bosons 
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We analyze the stability and decay of supercurrents of strongly interacting bosons on optical 
lattices. At the mean field level, the system undergoes an irreversible dynamic phase transition, 
whereby the current decays beyond a critical phase gradient that depends on the interaction strength. 
At commensurate filling the transition line smoothly interpolates between the classical modulational 
instability of weakly interacting bosons and the equilibrium Mott transition at zero current. Below 
the mean field instability, the current can decay due to quantum and thermal phase slips. We 
derive asymptotic expressions of the decay rate near the critical current. In a three dimensional 
optical lattice this leads to very weak broadening of the transition. In one and two dimensions the 
broadening leads to significant current decay well below the mean field critical current. We show 
that the temperature scale below which quantum phase slips dominate the decay of supercurrents, 
is easily within experimental reach. 



I. INTRODUCTION 

Many-body physics of strongly interacting ultra-cold 
atoms in optical lattices has been actively explored in the 
recent years. In particular, quantum effects of strongly 
interacting bosons, such as number squeezed states gen- 
eration^ and a quantum phase transition from the super- 
fluid to the Mott insulator—, have been observed in agree- 
ment with earlier theoretical predictions^. These devel- 
opments were followed by a broader theoretical analy- 
sis of phase diagrams of more complex systems includ- 
ing multi component bosongSiLSiSiH, Bose-Fermi mix- 
tures^^'^^'^^, and exotic states exhibiting topological or- 
j^gj.gi3,i4 gyg}j studies are motivated by issues that arise 
in a traditional condensed matter context. However, 
unique features of the cold atom systems also raise a com- 
pletely new set of questions. In particular, the ability to 
continuously vary interaction parameters, coupled to the 
near perfect isolation of these systems, open the way to 
address quantum dynamics far from equilibrium. 

In this context, there exists a purely dynamical phase 
transition of a condensate of weakly interacting bosons 
moving in an optical lattice. If the wave vector associ- 
ated with the condensate momentum exceeds a critical 
value, which is equal to one quarter of the reciprocal lat- 
tice constant for the square lattic o^^d^ , then the coherent 
motion of the condensate becomes unstable, resulting in 
the loss of superfluidity. Such a dynamical instability was 
observed experimentally^^ by measuring loss of coherence 
as a function of the condensate momentum. This tran- 
sition is of classical origin, in the sense that it is seen as 
an instability in the Gross-Pitaveskii equations of motion 
(GPE). Reltaed nonlinear dynamical phenomena such as 
self trapping and soliton formation have been studied in 
theory and experimentiSiiSiSfliSi. However, these studies 
focused on essentially classical systems, well described by 
GPE. Very little progress has been made in analyzing far 
from equilibrium behavior of systems where strong in- 
teractions and quantum fluctuations play an important 
role. 

In the present work we address this issue by focusing 



on a problem relevant to recent experiments: the fate 
of superfluid currents in optical lattices in the strongly 
interacting regime. We shall consider this issue via two 
questions, which will turn out to be closely related. First, 
what is the effect of quantum (as well as thermal) fluc- 
tuations on the dynamical instability of a moving con- 
densate? May the instabilitiy occur earlier, for example, 
in this case than the GPE prediction of pc — 7r/2? The 
second question is, how does the superfluid to Mott insu- 
lator transition take place when the condensate is moving 
in the lattice. This question may have direct relevance 
for ongoing experiments on the superfluid-insulator tran- 
sition. When the condensate is loaded into a magnetic 
trap it is hard to completely avoid center of mass oscilla- 
tions. In the absence of the optical lattice such a motion 
is frictionless and can persist for very long times. The 
same is true in the superfluid phase in the presence of 
the optical lattice as long as the center of mass momen- 
tum remains small and the interactions are weak. But 
what is the ultimate fate of this motion as the optical 
potential is increased and the system approaches the in- 
sulating regime? 

Effect of weak quantum fluctuations on the modula- 
tional instability in one dimensional traps was analyzed 
earlier numerically by one of ua^. It was shown that the 
quantum fluctuations smoothen the sharp classical tran- 
sition and lead to the current decay at smaller amplitude 
of the center of mass oscillations than predicted using the 
classical Gross-Pitaevskii (GP) equations alone. Similar 
numerical results were also reported in Ref. . Recent 
experiments confirmed strong damping of the center of 
mass oscillations in one dimensional condensates far from 
the classical modulational instability^^. 

In this paper we present a general theoretical frame- 
work of a superfluid-insulator transition in the current 
carrying state. Strictly speaking this is a true phase 
transition only at zero current. However, we find regimes 
where the broadening of the transition is small and even 
where a true discontinuity survives. In such cases a phase 
boundary is still well-defined. 

We shall show that at any nonzero current the transi- 
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tion is irreversible. If one starts from a condensate with 
nonzero current, increases the lattice strength past the 
transition point, then decreases it back to the original 
state, the current in the final state will have vanished. 
The energy contained in the initial motion of the conden- 
sate is transferred into thermal incoherent excitations. 
Thus the dynamical transition is of the first order type. 
The current carrying state is a metastable minimum of 
the classical (saddle-point) energy, and the transition oc- 
curs when the system escapes from this state. 

Throughout this work we employ the well known boson 
Hubbard model^^, described by the Hamiltonian: 

7i= ^ -J{a]ak + alaj) + ^^a]aj{a]aj ~ 1), (1) 

<]k> j 

where a] and aj are the boson creation and annihilation 
operators at the lattice site j, < jk > denotes pairs of 
nearest neighbors, J is the single particle hopping ampli- 
tude and U is the on-site repulsive interaction. Another 
implicit parameter is the average number of bosons per 
site N. If also a condition UN ^ J is fulfilled, then 
the boson Hubbard model can be mapped into the 0(2) 
quantum rotor model2^: 

n=Y^ -2JiVcos((/.fe - (/.,) + (2) 

<jk> j 

where (f)j and nj are the conjugate phase and the number 
of particles on the site j. The system described by ^ 
also undergoes a superfiuid insulator transition at JN ~ 
U and can support current in the superfiuid phase. In 
many situations, the quantum rotor model is easier to 
analyze analytically and we will frequently appeal to it. 
It is usually well justified in the case ^ 1. Indeed, in 
this limit it is possible to have simultaneously UN ^ J 
and JN either smaller or larger than U. So at large N 
the mapping from the boson Hubbard to the quantum 
rotor model can be justified in both the superfiuid and 
insulating phases. 

Our paper is organized as follows. In sec. ^ we give an 
overview of our main results and present a physical pic- 
ture. In sec, mil we derive the mean field phase diagram 
separating stable and unstable regimes of current flow. 
Sec. IIVI focuses on the current decay mechanisms due to 
quantum and thermal fluctuations. In particular, we ob- 
tain leading asymptotic contributions to both quantum 
and thermal decay rates near the mean field instability. 
Then in sec.^we consider dynamics of the current decay 
and discuss the effects of a parabolic confining potential. 
Sec. IVII addresses the loss of coherence in a condensate 
following super-current decay. In sec. lVIII we present the 
results of exact simulations in small systems and discuss 
them in the context of our theoretical analysis. Finally, 
in sec. IVIlH we summarize the results and discuss exper- 
imental implications. A shorter discussion of the results 
presented here can be found in Ref. (26| . 



II. PHYSICAL PICTURE AND OVERVIEW OF 
THE RESULTS. 

The existence of a critical velocity of a condensate in 
an optical lattice was predicte(ii^*i& and observed exper- 
imentally 2Li^ in the regime of weak quantum fiuctua- 
tions {JN ^ U). In this case one can solve the Gross- 
Pitaevskii equations of motion and find that the conden- 
sate becomes unstable when the phase difference between 
adjacent sites becomes larger than tt/2. There is a simple 
way to understand this instability by considering how the 
superfiuid current flowing through the system changes 
with the condensate momentum. In a continuum system 
the current is just 

J = PsP, (3) 

where ps denotes the superfiuid density and p the mo- 
mentum or phase gradient such that (j){x) — px. In a 
discrete system described by a tight-binding model, the 
above expression is modified to 

J ^ Ps sin p. (4) 

More generally sinp is replaced by a different function 
of p, with the same elementary period. At small cur- 
rents we recover the continuum limit from Eq. In 
the Gross-Pitaevskii regime, the superfiuid density does 
not depend on momentum. Therefore the maximal cur- 
rent occurs at p — 7r/2, precisely where the condensate 
motion becomes unstable. As we argue below, this is not 
a coincidence. Consider a perturbation, where the state 
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FIG. 1: Schematic representation of a perturbation to a state 
with a uniform phase gradient. Dots represent phases on dif- 
ferent sites for uniform and perturbed systems. The lines are 
guides to an eye. 

with a uniform phase gradient p is split into two equal 
domains with slightly higher and slightly lower momenta 
pzt Sp (see Fig. At small Sp we can expand the en- 
ergy difference between the perturbed and unperturbed 
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configurations in powers of 5p. The linear term in the 
expansion vanishes because the contributions to the en- 
ergy from the left and the right domains exactly cancel 
each other and we are left only with the quadratic term: 

IdE^ IdE, ol'^^^/cN2 l'^^^.. ^2 

" 2 2 d^(-'^) + 24 - 2 V ^'"^ ' 

(5) 

Noting that the superfluid current is formally defined as 
the derivative of the energy with respect to phase gradi- 
ent: J{p) — dE/dp we find 

SEoi^iSpf. (6) 
dp 

Thus if the current is an increasing function of the mo- 
mentum, then small deviations from the uniform state 
increase its energy. On the other hand if dj /dp < 
then the fluctuation just considered reduces the energy 
of the system. In this case it is obvious there is a manifold 
of resonant configurations obtained by a smooth contin- 
uum transformation from a uniform state. For example 
one can take the state in Fig. ^ plus a weak positive en- 
ergy fluctuation. So there is no energy barrier protecting 
a uniform state from fragmentation. Let us note that 
this argument shows that the positive slope of the cur- 
rent with respect to p is a necessary condition for the 
stability of the uniform state. It does not exclude, how- 
ever, that the current can decay even if this condition is 
satisfied. In sec. II VI we will see this indeed occurs due to 
quantum or thermal phase slipsSLS^. 

Consider now the strongly interacting regime. Suppose 
that we deal with a uniform system at commensurate (i.e. 
integer) filling. It is well known that such systems un- 
dergo a superfluid-insulator transition^ at zero tempera- 
ture. This transition is driven by quantum fluctuations 
which increase with the interaction strength. As p in- 
creases the effective hopping amplitude in the direction 
parallel to the current decreases as Jeff J cosp result- 
ing in reduction of the sound velocity^®. Alternatively 
reduction of Jg/ / can be viewed as increase of the single- 
particle effective mass with quasi-momentum. This im- 
mediately follows from a single-particle band structure. 
As a result, quantum fluctuations, which are determined 
by the ratio (U/JeffN) |3^, become stronger with p, im- 
plying concomitant increase of quantum depletion of the 
superfluid density. Therefore the equation 10} should be 
rewritten as: 

J = p(p) sinp, (7) 

where p{p) is a monotonically decreasing function of the 
momentum. Thus the product reaches a maximum at 
some p* < Tr/2. In the Gross-Pitacvskii regime JN ^ U, 
the dependence of p{p) on p is very weak and we find 
that p* ~ tt/2. On the other hand, in the vicinity of the 
superfluid-insulator transition p{p) is both very small and 
very sensitive to variations of Jef / . Thus we expect that 
in this case p* will be close to zero. 



In sec mil we give a formal derivation of the critical 
momentum at which the condensate motion becomes un- 
stable. Using the time dependent Gutzwiller approxi- 
mation, we plot the critical momentum as a function of 
the interaction strength. This phase boundary interpo- 
lates between the usual dynamical instability occuring at 
p = ■k/2 for small interactions and the vanishing critical 
momentum at the equilibrium superfluid-insulator tran- 
sition (see the top graph in Fig. 0}. 

The situation is different in the non-commensurate 
case . No matter how strong the interaction strength, 
it can not localize the excess particles (or holes) moving 
on top of the Mott background. This excess density al- 
ways remains superfluid, independent of Jeff and thus 
also of p. Together with Eq. {T)), this implies that at 
strong interactions the instability occurs at p = 7r/2. On 
the other hand for sufhciently weak interactions, where 
the number fluctuation per site SN >> 1, there is no dis- 
tinction between integer and non integer density. There- 
fore for U not too large the critical momentum should 
decrease with U from 7r/2. Indeed we find, using the 
time dependant Gutzwiller approximation, that for the 
incommensurate case Pc reaches a minimum at some fi- 
nite interaction strength and saturates on 7r/2 for both 
very weak and very strong interactions, (see the bottom 
graph in Fig. Q}- 

In Sec, nil 51 we develop an analytical approach, which 
describes the behavior of the critical momentum pc ver- 
sus interaction in the vicinity of the zero-current SF-IN 
transition. We show that in the commensurate case Pc 
vanishes as 1/^, where f is a coherence length, which di- 
verges at the transition. At non-integer filling we confirm 
the non monotonic behavior of the critical momentum. 

In practice the system always includes a global har- 
monic confinement, which leads to a non uniform density 
distribution. In this case we find that the instability first 
occurs at the the edges of the condensate where = 1 
regardless of the peak density A'o in the middle of the 
trap (see Sec. 01. There is a difference between large and 
small Nq, which manifests in the dynamics of the current 
decay. For Nq « 1 (as well as in uniform systems with ar- 
bitrary filling) we find that near the instability the decay 
is underdamped, i.e. the instability rapidly grows in time 
destroying the current state. On the other hand if Nq is 
large, the momentum oscillations decay gradually after 
the edges become unstable. There is another important 
difference between the uniform and parabolically trapped 
condensates. In a uniform lattice there are only two en- 
ergy scales, set by U and J and their ratio completely 
determines the behavior of the system. With harmonic 
confinement on the other hand, due to the presence of 
another energy scale one should take into consideration 
whether [/ or J or both are being changed in the exper- 
iment (see discussion in Sec. FVland Appendix 0. 

So far we concentrated on the results of the mean field 
dynamics at zero temperature, where the time evolution 
is simply described by classical equations of motion. In 
Sec. lIVI we go beyond the mean field dynamics to analyze 
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the effect of quantum and thermal fluctuations. These 
act to generate phase-slips, which induce current decay 
even prior to the classical instability. A phase config- 
uration for a particular phase slip is shown in Fig. |21 
Basically a phase-slip corresponds to generation of a large 
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FIG. 2: Schematic representation of a phase-slip. Notations 
are the same as in Fig. Q 

phase difference on a particular link (or in the vicinity of 
this link) and simultaneous reduction of the phase gradi- 
ent in the rest of the chain. Because the energy is a pe- 
riodic function of phase differences, by generating phase 
slips the system reduces its superfluid current. We cal- 
culate the leading exponents of the decay rates in the 
Gross-Pitaevskii regime of relatively weak interactions 
and in the quantum critical regime close to the SF-IN 
transition. We find, that broadening of the mean field 
transition is strongest in the one-dimensional case. In 
particular, deep in the superfluid regime, the phase-slip 
tunneling rate at p — > 7r/2 scales as 



r(xexp(^-5,^:i^(7r/2-p)^j , (8) 

where Sd is just a number that depends on the number of 
space dimensions (see Eqs. We obtain similar 

results for the thermal decay. Although below 7r/2 the 
decay in 3D is clearly much weaker than in ID, there is 
no qualitative difference between various dimensions. As 
long as the ratio JN/U remains large, the tunneling of 
phase-slips is exponentially suppressed. 

Next we derive analytical expressions for the exponents 
characterizing the decay rate in the vicinity of the equi- 
librium SF-IN transition. We find again that fluctuation 
induced decay is stronger in lower dimensions. However, 
there is no small parameter like U/JN controlling the 
mean field results. We show that in one dimension the 
exponent always remains of the order of one, implying 
significant broadening of the mean field transition. In 
three dimensions, by contrast, we find that the quantum 



decay rate does not vanish at the mean field instability, 
but rather exhibits a discontinuous jump. In this sense, 
the three dimensional system undergoes a sharp dynam- 
ical transition at zero temperature. 

The physical picture of current decay due to generation 
of phase slips is similar to the situation in superconduc- 
tors. In particular, deviation of the critical current from 
the mean field result was observed in thin superconduct- 
ing wires^^ and explained theoreticallji^Si^^ . The mecha- 
nism responsible for reduction of the critical current was 
identified as creation of phase slips due to thermal fluctu- 
ations. The question of observing current decay in super- 
conductors due to quantum tunneling is still debated (see 
Ref. for recent developments). We will show that in 
the systems under consideration here, current decay due 
to quantum effects is easily within experimentally reach. 
For a one dimensional system, for example, the character- 
istic temperature below which quantum decay dominates 
in the GP regime (far from the Mott transition) is of the 
order of the Josephson energy (T* « VUJN). This is 
much higher than typical temperatures in optical lattice 
experiments. At strong interactions, in the vicinity of 
the Mott transition, broadening of the mean field tran- 
sition due to quantum phase slips is always large in one 
and two dimensions. Only in the three dimensional case, 
where quantum tunneling is suppressed, thermal fluctua- 
tions can be responsible for current decay below the mean 
field transition. Some additional details on the relation 
between current decay in superconductors and in optical 
lattices can be found in Ref. js^ . 

Let us now briefly mention some interesting experimen- 
tal implications of our results. We envision the following 
experimental scheme. Start with a superfluid state far 
from the Mott insulator. Then either boost the conden- 
sate to some non zero momentumii, or induce a center 
of mass oscillation in the trap24. Now, slowly increase 
the interaction parameter up to a specified point. This 
can be done in the usual way by increasing the lattice in- 
tensity, or alternatively by decreasing the detuning from 
a Feshbach resonance. Finally decrease the interactions 
back to their original value. If the dynamical phase tran- 
sition is sharp then as long as the system does not cross 
the transition boundary (Fig. ^) within this cycle, the 
process should be completely reversible. In particular, 
the initial current (or center of mass oscillation), as well 
as phase coherence, should be fully recovered. At the 
same time, if the system does pass through the transi- 
tion, the current will be lost irreversibly and the system 
will heat and partially lose its coherence, compared to 
the original state. 

One of our main results, is that in a three dimen- 
sional optical lattice this mean field dynamical transi- 
tion is sharp, and it essentially survives the effect of 
fiuctuations. Such experiments can thus map the non- 
equilibrium phase diagram shown in (Fig. and directly 
trace the connection between the classical modulational 
instability {pc = 7r/2) and the the equilibrium Mott tran- 
sition. In fact, due to it's discontinuous nature, the dy- 
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namical transition point is much easier to detect. This 
suggests an accurate method to determine the position 
of the equihbrium Mott transition by extrapolating the 
dynamical transition line to zero momentum. 

One comment is in order concerning heating and loss 
of coherence in the final state. In Sec. I VII we show that in 
three dimensions this loss of coherence is only significant 
at very large currents {p ~ 7r/2). In one dimension (and 
to a lesser extent in two dimensions) , the loss of coherence 
due to irreversible heating depends on the system size or 
experimental resolution and may thus be large even at 
small currents. 



III. MEAN FIELD DYNAMICS AND CRITICAL 
CURRENTS IN OPTICAL LATTICES. 

A. Gross-Pitaveskii regime 

Despite its simple form, the Bose Hubbard model is 
not integrable in any spatial dimension^^*^ and can not 
be solved completely. Nevertheless, there are some limits 
where one can make considerable progress in understand- 
ing its static and dynamic properties. In particular, one 
can easily address the regime of weak quantum fluctu- 
ations, which is the case when JN ^ U—. Then one 
can use discrete Gross-Pitaveskii equations^. For the 
Hamiltonian ^ these are given by: 



. dipj 



(9) 



where the classical fields ipj and ■0^ correspond to the 

expectation values of aj and respectively; the set O 
contains the nearest neighbors of site j . In the quantum 
rotor limit UN 3> J the number fluctuations are weak 
and can be integrated out leaving us with equations of 
motion for the phase 4)j = argi/^j only: 



-2UJN sin( 
feeo 



(10) 



Both equations © and H10|l can support stationary cur- 
rent carrying states ipj oc exp(zpj). A simple linear sta- 
bility analysis shows that these states are stable towards 
small perturbations for p < tt/2 and become unstable 
otherwis o ^ "'M^ . Thus, 7r/2 is the critical phase twist above 
which a uniform superfluid state breaks down. The tran- 
sition from the superfluid to the insulating state asso- 
ciated with this instability is known as the classical lo- 
calization transition. It was recently observed experi- 
mentally^^. In the presence of quantum fluctuations, the 
current can decay even for p < tt/2 via quantum tun- 
neling. Clearly these fluctuations should be increasingly 
important as the system approaches the Mott phase. The 
same is true for decay due to thermal fluctuations as one 
increases the temperature. In the next section we will 
address this issue in detail. 



B. Critical current in the vicinity of the SF-IN 
transition 

The Gross-Pitaveskii description of the dynamics 
breaks down at strong interactions. Moreover when 
JN ^ U the bosonic system at commensurate filling 
{N is integer) undergoes the Mott insulator transition 
entirely driven by quantum fluctuations^iS. In the uni- 
form system with a fixed density this transition lies in the 
universality class of the xy model in d -I- 1 dimensions'^, 
the properties of which were well studied long ago2^. So 
there is also a hope to get insights to some nonequilib- 
rium properties of the interacting bosons in the vicinity 
of the phase transition. The latter, as any generic sec- 
ond order phase transition, is characterized by a diverg- 
ing correlation length ^ |25|], which sets the length scale 
for all low-energy universal properties of the system. In 
particular, close to the critical point the low-energy long- 
wavelength fluctuations can be described by relativistic 
(z = 1) effective field theorjiS^. In terms of dynamics this 
implies that the classical equations of motion also take 
explicitly relativistic invariant form^^: 



(11) 



where ip is the superfluid order parameter, r tunes the 
system across the SF-IN transition: r > corresponds 
to the superfluid phase and r < does to the insulator. 
Here we rescaled the units of coordinates and time by 
a constant of the order of one (see Appendix ^ for the 
details). The correlation length ^ is related to the tuning 
parameter r as ^ cx 1/ ^/\r\■ We point out that Eq. 1)11(1 
is very reminiscent of the conventional continuum Gross- 
Pitaeskii equation with the only difference that there is a 
second (as opposed to first) order time derivative in the 
LHS. This equation has a conserved charge: 



(12) 



which is proportional to the deviation of the particle den- 
sity from the integer filling; the constant A is explicitly 
given in the Appendix ^ Thus the stationary solutions 
correspond to the commensurate transition. In the non- 
commensurate regime there is no phase transition, but 
one can still use Eq. 1)11(1 if the deviation from the inte- 
ger filling is small. 



1. Commensurate case 

Let us analyze the fate of the current-carrying case if 
the mean number of bosons per site is integer. Equa- 
tion (|11() supports stationary states of the form: 



(13) 



where z denotes d — 1-dimcnsional space of transverse 
coordinates. Without any further analysis it is obvious 
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that current states become unstable at p > \fr ^ Xj^^, 
i.e. the critical momentum vanishes at the superfluid- 
insulator transition point. To be more precise, we can 
evaluate the spectrum of small fluctuations of Eq. 
around the stationary solution (|13|l : 



u;2(q) 



V 



q2 ± V(r-p2)2+4p2^2^4^ (14) 



where q is the wavevector characterizing the fluctuations 
around the state (|13|l . In the long wavelength limit |q| 
the expression above yields simplified frequencies for 
the amplitude and the phase modes: 



cc>2(q)«2(r-p2) + 



-2 2 



'^2(q) 



3p2 



2 



(15) 
(16) 



The first (amplitude) branch is always gapped unless 
> r and therefore is stable against small fluctuations. 
On the contrary the second, phase mode, becomes unsta- 
ble at p > pc = We would like to stress that the 
relativistic nature of excitations is crucial to get this in- 
stability. Despite being continuum equations (|ll|l rely on 
the presence of the underlying lattice, which breaks the 
translational invariance. Otherwise the equations of mo- 
tion would be Gallilean invariant and no critical current 
would exist. 

From the analysis above we see that close to the 
superfluid-insulator phase transition current states be- 
come unstable at momenta inversely proportional to the 
correlation length of the condensate. As one goes deeper 
into the superfluid regime the correlation length de- 
creases saturating at one (in the lattice units) and we 
come back to the Gross-Pitavcskii result of instability 
occurring atp = 7r/2~l. 



2. Incommensurate case 

It is also interesting to consider the stability of the 
current states at the noncommensurate filling. In this 
case the system remains superfluid at arbitrarily strong 
interactions^. If the interactions are weak the system 
is in the Gross-Pitaveskii regime and the flUing is not 
important. In this case we expect a usual modulational 
instability at p «:! 7r/2. At the same time, when the inter- 
action strength becomes very large, we can think about 
excess particles as hard-core bosons moving on top of the 
Mott insulator. But in turn, this would be equivalent to 
a spin one half system. At the mean fleld level we can 
use again spin- wave theory to see that pc is exactly equal 
to 7r/2. This suggests that pc should have a minimum at 
the intermediate interaction strength. 

The solution of Eq. I|ll(l corresponding to the noncom- 
mensurate filling factor can be written as: 



ip{x, z, t) — pe 



(17) 



where p = y/r + p? — p'^ and p is related to the deviation 
from integer filling Sn: 



5n = App^ 



(18) 



As in the commensurate case, in the long wavelength 
limit there are two branches describing a gapped ampli- 
tude mode and gapless phase fluctuations. The disper- 
sion of the latter for q parallel to the x-axis reads 



uj{q) 



2mP 

2^2 + p2 - ■ 2p- + p 



■q 



+ . /, , V3p^-2r|g|. (19) 



From this we observe that the current state first becomes 
unstable when Bp^ = 2r. Together with l|18|l this gives 
the critical momentum 



Pc 



9 5v? 



(20) 



This result reduces to the commensurate limit for 5n = 0. 
However, for any nonzero 5n it shows that pc reaches 
the minimum value pj oc Sn^^^ at rg oc 8r?l^ and then 
diverges as r ^ 0. While the divergence is the spurious 
result of the continuum approximation; it should be cut 
of by the lattice at p w 1, the existence of the minimum 
agrees with the simple general argument given above. 



C. Gutzwiller approximation 

Having derived the conditions for the stability of the 
current-carrying condensate in a lattice in the two ex- 
treme limits, one can try to find a unifying approach, 
which interpolates between them. A natural choice is 
the Gutzwiller approximation. This is a time-dependent 
generalization of the standard mean field theory, where 
the wavefunction is assumed to be factorizable: 



\G\ 



n 



(21) 



Here j denotes a site index and n site occupation. The 
ansatz ()21|l supplemented by self-consistency conditions 
leads to equations of motion: 



-Jz{Vnfj^n-i'fpj + Vn-I- 1/j-n+i?/'*), (22) 



where 



(23) 



ieo 



O is the set of nearest neighbors to j and z is the co- 
ordination number (z = 2c? for a hypercubic lattice). In 
practice the sum in H21|l is limited to a finite number 
of states, on each site, so that only a flnite number of 
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equations need to be solved. We checked that in our nu- 
merical simulations we take sufficient number of states 
so that the results are insensitive to the truncation. In 
particular, for N = 1 we compared simulations for the 
spectrum truncated at 5 and 10 states per site and the re- 
sults were practically indistinguishable. The Gutzwiller 
approximation can be justified at high dimensions, where 
the coordination number z becomes large. In this sense 
it is reminiscent to the dynamical mean-field theorji^. 
In reality, it is necessary to calculate first quantum cor- 
rections to the evolution governed by Eq. (|22|l to see the 
validity of the Gutzwiller result at a given dimensionality. 
We will postpone such an analysis until the next section. 

To find numerically the position of the dynamical in- 
stability corresponding to Eq. (|22|1 we can carry out one 
of the following procedures, (i) Starting from the non- 
interacting state {U = 0), where the Gutzwiller ansatz 
becomes exact, and a given phase gradient p we adia- 
batically increase U. Observing either the current or the 
condensate fraction (which we define as the population 
of the state with the momentum p) we can identify the 
critical interaction at which the motion becomes unstable 
(see Fig.O. (ii) Alternatively we can find numerically the 




I -I ' 1 ' 1 ' 1 ' 1— 

0.0 0.2 0.4 0.6 0.8 

Interaction U/U 

c 



FIG. 3: Condensate fraction as a function of scaled interac- 
tion for a two-dimensional condensate with filling N = 1 eval- 
uated within Gutzwiller approximation on a square lattice of 
size 80 X 2 |4lll . The initial momentum is p = 7r/5. The other 
parameters are J = 1 (so that the interaction strength corre- 
sponding to the transition to the Mott state is Uc ~ 23.2) and 
the interaction increases in time as U = O.OAt. The conden- 
sate motion clearly becomes unstable at certain interaction 
{U/Uc ~ 0.57), which marks the dynamical transition. 



mean field ground state corresponding to given U and J 
and adiabatically increase a gauge potential so that H23|) 



modifies to 

V-.M^ i(^(G|aj+i.i|G)e''^(*) 

+ (G|a,_i,,|G)e-^'^W-t- ^ (Gla.-HG)]. (24) 

I'GO' ^ 

Here we explicitly introduced indices along the current 
j and in the transverse direction 1; O' is a subset of O, 
which excludes the sites {j ± 1, 1}. If the system is stable 
then the condensate remains at the momentum p — 
in a moving lattice, which is of course equivalent to a 
moving condensate with p = (j){t) in a stationary lattice. 
Once the motion becomes unstable the distribution at 
p — rapidly drops. The second approach is favorable, 
because it does not require quantization of the momen- 
tum in units of 27r/M in finite systems of longitudinal 
size M, where the actual calculations are performed. 

Identifying the point of dynamical instability as de- 
scribed above for different interaction strengths we can 
construct a mean field phase diagram separating stable 
and unstable condensate motion for both integer and 
noninteger filling factors (see Fig. This phase dia- 
gram is in complete agreement with what we expected 
from the analysis given in the previous subsections. Thus 
at small interactions the critical momentum approaches 
7r/2 for any filling or dimensionality of the lattice. At in- 
teger filling the critical momentum vanishes at the point 
of commensurate superfluid-insulator transition. In the 
incommensurate case the critical momentum first goes 
down with the interaction strength and then increases 
back to 7r/2 in the strongly interacting regime. 



IV. BEYOND MEAN FIELD THEORY. 
CURRENT DECAY DUE TO FLUCTUATIONS. 

The analysis given in the previous section is valid only 
at the mean-field level. Quantum and thermal fluctua- 
tions can destroy the current by either phase slip tunnel- 
ing or thermal activation. The main goal of this section 
is to derive leading contributions to the decay rate and 
check the validity of the mean field results. To simplify 
the analysis we will concentrate on the two tractable lim- 
its: the Gross-Pitaveskii regime describing the system 
deep in the superfluid phase and the Ginzburg-Landau 
regime, which is valid in the vicinity of the superfluid in- 
sulator transition, where the correlation length becomes 
large compared to the lattice constant. Also for simplic- 
ity we assume that the filling is large and integer so that 
one can use the quantum-rotor model. 
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FIG. 4: Mean field phase diagram separating stable and un- 
stable motion of condensate. The vertical axis is the conden- 
sate momentum in the inverse lattice units and the horizontal 
axis is the normalized interaction. The top graph shows the 
result for integer filling = 1 at different spatial dimensions. 
The bottom graph describes a two-dimensional lattice with 
different filling factors. 



A. Gross-Pitaveskii regime. 

1. Current decay due to quantum tunneling 

As we argued above, the current state with p < 7r/2 is 
stable with respect to small fluctuations. However, this 
state does not correspond to the energy minimum, which 
has no current. So we conclude that such a state must 
be metastable. Contrary to a uniform system, where 
the momentum can be gauged away via the Gallilean 
transformation, in a lattice there is a preferred reference 
frame. This immediately implies that such a metastable 
state should be able to spontaneously decay because of 
quantum tunneling. 

In the leading order in U/ JN, which plays the role 



of the effective Planck's constant for this problem'*^, the 
tunnehng rate (F) corresponds to the action {S) of the 
bounce solution (instanton) of the classical equations of 
motion in the inverted potential^: 



F cx e" 



(25) 



We will not attempt to compute the prefactor here and 
will concentrate only on S. We just point out that the 
prefactor scales with the system size and in the thermo- 
dynamic limit the tunneling rate per unit volume is size 
independent. Explicitly the action reads as: 



1 



2U \ dr 



2JN cos(0j4+i 



2JNcos{p + - 



(26) 



where as before j denotes the coordinate along the cur- 
rent direction and 1 corresponds to the transverse co- 
ordinates. In H2t)|) we redefined phases compared to 
H10|l : (fij 4>j — pxj so that ipj = corresponds to 
a metastable current state and the new phases are the 
fluctuations around this state. It is convenient to rede- 
fine the imaginary time r, measuring it in Josephson time 
units: t ^ r/ VuWj. Then it is easy to see that 



S=^l—s, 



(27) 



where 



1 fdcf) 



2 cos((?!)ja+i - ipj^) 



2 V dr 

2coa{p + ipj+i^i - 



(28) 



The desired instanton trajectory is the solution of the 
Euler-Lagrange equations, which can be obtained ex- 
tremizing the action with respect to the phases 0j i and 
subject to a boundary conditions: = at r = ±oo 
and at |1|, \j\ — oo. 

Before proceeding with a general analysis in higher spa- 
tial dimensions let us consider the case d = I first. In 
this work we are interested in the decay close to the crit- 
ical current pc = tt/2. Clearly as p ^ Pc the effective 
tunneling barrier becomes weaker and weaker and hence 
we can expand (|28|) in powers of : 



dr 



■cos(p) {(f) j+i - 0j)" 



.(29) 



This expansion is similar to that used in the analysis of 
thermally mediated decomposition near spinodal pointM. 
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In the action above we used that sin p « sin pc = I. Now 
we can do another rescahng: 



— COs(p) (f). 



y/cOs{pj' 

which simphfies the action even further: 



(cosp) 



5/25 



(30) 



(31) 



where 



df 



+ 



Pj+i 



(32) 

Note that s is just a number of the order of one so that 
(|31|l completely determines the momentum dependence 
of the action. The scaling (|30|l guarantees that the orig- 
inal phases 0j remain small for the instanton trajectory 
and H29|l is indeed a correct asymptotical form of H28|l . 

Since the action s contains no small parameters, the 
bounce solution should be localized within a few sites. 
Without loss of generality we can assume that the maxi- 
mum phase difference develops between the sites labeled 
by i = and j = 1. For the rest of the system the 
phase gradients are relatively small so we can neglect 
cubic terms. Then those degrees of freedom can be inte- 
grated out: 



An 



(33) 



where a^r) = (r) — (/)o (r) and a{uj) is its Fourier trans- 
form; 



(34) 



Substituting into HT^ we find: 



dT 



1 / da{i 



4 V 

i 47r' ^ ^' 1,^1 



1 



+ a^(T) - - a'^(T) 



V8 



(35) 



Clearly, if we ignore the last term in the expression 
(|35|l we get the action of a single particle moving in 
a metastable potential. The last term represents a 
dissipative-like contribution coming from the rest of the 
chain. If we ignore this term then the solution extremiz- 
ing is: 



a{T) 



cosh^ t/\/2 



(36) 



This yields s = 24/5 = 4.8. In a general case with dissi- 
pation we can try a variational ansatz solution: 



q;(t) 



A 



cosh rr 



(37) 



A simple numerical analysis gives 

r « 0.64, A « 3.31, S « 7.11. 



(38) 



So the action is about a factor of 1.5 larger than without 
the bath degrees of freedom. Using (|^ and one 
can verify the consistency of the harmonic approximation 
used for the sites other than "1" and "0". In particular, 
it is straightforward to get: 



102(0) -01 (0)1 
«(0) 



0.326, 



(39) 



which is relatively small. The difference between phases 
in further nearest neighbor sites is even less. Therefore 
for them the harmonic approximation is justified even 
better. 

Instead of harmonic treatment of the phases other than 
and 02 we can exactly take into account the four near- 
est neighbor sites and ignore the rest of the chain. Then 
the instanton solution is parametrized by the two inde- 
pendent angles a and f3: 



PO = 



= a/2, 4)1 = -0-2 = a/2 + /3. 



(40) 



Substituting this into the action and solving the corre- 
sponding Euler-Lagrange equations one can show that in 
this case s ~ 6.1, which is about a factor of 1.26 larger 
than the result with /3 = 0. This number is the exact 
lower bound for the action s since the other degrees of 
freedom can only increase the action. We will not fur- 
ther attempt to improve the accuracy of s noting only 
that the variational result s « 7.1 should be very close 
to the exact value. 

We can straightforwardly generalize the one- 
dimensional results to higher spatial dimensions. 
In particular, using the same arguments as before close 
to the critical current we can expand the action (|28|l up 
to the cubic order in (pjy. 



\[^] + cos(p) (0j+ia - 0,- 1)^ 



I'eo' 



i^3-\f - 3(0^+1,1 - 0j,i)^ 



(41) 



Note that at p — > 7r/2 only longitudinal modes become 
soft, acquiring a prefactor cosp in front of the quadratic 
term in the action. This implies that the transverse width 
of the instanton should be much larger than its longitudi- 
nal size and we can safely use the continuum approxima- 
tion for the phases along the transverse directions. Then 
instead of H41|) we derive: 



drd 



dT 



dx 



+ COs{p) {(t)j + i ~ <l>jf - -( 



(42) 
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In this equation x denotes transverse coordinates which 
reside in a d — 1 dimensional space. As in the one dimen- 
sional case we can rescale 



4> = COs(p) (j), T 



T 

^cos(p) ^cosp 
In this way the action (|42|l becomes: 

where 



(43) 



(44) 



c 




+ 



1.7 



Vj+l, 



(45) 



Here C = [t, x) is a d-dimensional space-time coordinate. 
As before Sd is just a number, which depends only on di- 
mensionality. The result H44I) is quite remarkable. First 
of all it shows that the action in higher dimensions van- 
ishes much more slowly near the critical current. From 
the scaling (|43|l it is obvious that the characteristic trans- 
verse dimension of the instanton scales as 1/VPc — P 3* 1, 
justifying the continuum approximation. Above d = 6 
the tunneling action would experience a discontinuous 
jump at p — Pci however since we deal with d < 3 this 
is not relevant. Now let us try to estimate Sd- We again 
proceed in the same spirit as in the one-dimensional case. 
In the first approximation we consider only a single phase 
slip (pi = a/2, (j)2 — —a/2 and treat the motion of other 
phases in the harmonic approximation. The correspond- 
ing dimensionless action reads 



Sd = 

1 

+ 2 



d-1 

2 — 



1 3 

3" 



(2^' 



.(k)l^ 



(46) 



where a(k) is the Fourier image of a(x). If we ignore the 
last dissipative term in l|46|) . then our action is identical 
to that considered in Ref. ,43|] for the decay of a false 
vacuum. In that work it was argued that the bounce 
solution is spherically symmetric and it satisfies the fol- 
lowing equations of motion: 



1 J__d_ f.d^\ 



2a 



a 



(47) 



with the boundary conditions: a(oo) = and a'(0) = 0. 
These equations can be easily solved numerically and the 
result is: 



h ~ 21.92, S3 ~ 87.32. 



(48) 



So it is clear that at higher dimensions not only the ex- 
ponent of (pc — p) in the instanton action gets lower but 



also the numerical prefactor gets larger. Now we can 
find the correction to the action due to the bath degrees 
of freedom coming from the last term in 1)46(1 . Instead 
of using the variational approach as we did in the one 
dimensional case, we will use the exact solution of H47|l 
to evaluate the contribution of the bath term in the ac- 
tion. This will be the exact upper bound of the action 
Sd- Direct evaluation of (|46|l gives 

4.8 < si < 8.0, 21.9 < S2 < 27.6, 87.3 < h < 97.5. 

(49) 

Obviously the local approximation = for j ^ 0^1 
works better and better at higher dimensions implying 
that the effective size of the instanton in the longitudinal 
direction (along the current) decreases with the dimen- 
sionality of the space. 

Because the decay rate is strongly (exponentially) de- 
pendent on momentum and coupling constants we can 
approximately define the stable phase at which the tun- 
neling action is larger than some number, say S > 3 
and unstable phase, when there is no exponential sup- 
pression of the tunneling of phase-slips, say S < 1. The 
region in between will denote the crossover between the 
stable and unstable regimes. In this way we can define 
a crossover phase diagram (Fig. (SJ . We see that except 



0.50 



o 0.45 4, 
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FIG. 5: Large A'' zero-temperature phase diagram for the 
nonequilibrium superfluid-insulator transition. The dashed 
line represents mean field transition at p = n/2. The shaded 
regions correspond to the tunneling action satisfying 1 < S < 
3, obtained within the discrete phase model in diflterent spatial 
dimensions. Below the shaded regions the tunneling action is 
large and the current decay is slow, so the superfluid state 
is stable for long time scales. Above the shaded regions the 
current decays very fast and the superfluid state is unstable. 

for U/JN <C 1, there is a very strong broadening of the 
classical transition in one dimension. On the contrary 
in the three-dimensional case effects of quantum fluctua- 
tions are relatively weak and the crossover is very sharp. 
We would like to point out that the derivation given here 
is valid only deep in the superfluid regime U/JN <C 8d. 
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Close to the critical point it is necessary to use the coarse- 
grained description which we discuss later. 

To summarize this subsection we write explicit expres- 
sions for the tunneling action in the phase model in three 
spatial dimensions: 



~ 25 



(50) 
(51) 
(52) 



2. Thermally activated current decay. 

Now let us turn our attention to broadening of the 
mean field transition due to thermal fluctuations. A 
general formalism for finding the decay rate was devel- 
oped by Langer^^. It was later successfully applied to 
quasi-one-dimensional superconductors^Si^ and to three 
dimensional superfluids at small currents^. 

Before proceeding with this general method let us 
point out an essential difference between conventional 
condensed matter systems and the cold atoms systems 
addressed in this work. In the former, it is the environ- 
ment, which introduces thermal noise and dissipation**^ 
to the system leading to diffusion in energy space and 
eventually thermal activation. In cold atom systems by 
contrast, the temperature is introduced into the system 
during the preparation of the condensate. I.e., the initial 
state of the condensate is described by a thermal den- 
sity matrix rather than a pure wave function. Later the 
system is essentially isolated from the environment and 
evolves according to the Hamiltonian equations of mo- 
tion. We point out that the formalism of Ref. was 
based on very general assumptions, and therefore should 
be independent of the details of the thermal fluctuations. 
Nevertheless, certain issues arise that require special at- 
tention. Consider for example a superfluid flowing in a 
container whose walls act as the thermal bath. The wall 
as well as the thermal fluctuations arising from it set a 
preferred reference frame, breaking the galilean invari- 
ance of the superfluid and thus allowing for the current 
decay. An isolated superfluid, on the other hand, even 
if prepared at flnite temperature, is galilean invariant. 
Thus current in such a superfluid cannot decay unless 
there is an external potential, such as a lattice, that sets 
a preferred reference frame. Because we are interested in 
the current decay in the limit where the lattice is strong 
and BHM is applicable, this subtlety is irrelevant for our 
consequent analysis. The effect of breaking of galilean 
invariance by weak external potential on thermally acti- 
vated current decay in one dimensional superconductors 
was recently studied in Ref. jiS ]. 

To simplify the analysis we assume that the temper- 
ature is small so that there is no difference between the 



energy and the free energy. Indeed, the difference be- 
tween the two amounts to the product TS, where S is 
the entropy of the system. At small temperatures the 
latter vanishes in the superfluid as T'^ so that the prod- 
uct TS goes to zero at least as at T — + and is always 
negligible compared to the activation energy, which does 
not depend on temperature. 

As in the previous subsection, we first c onsider the 
one dimensional case. Following Refs. |3ll45l | we find the 
stationary solutions of the classical equations of motion: 



1 



2 dr^ 



sin(p + (/)j+i — 0j) +sm{—p + (f>j^i — (f>j). (53) 



Clearly = describes a metastable state carrying 
the current 2J7Vsinp. We note again that phases in lfS5|l 
are the deviations from the metastable current state. The 
other saddle point solution separating the states with dif- 
ferent currents is 



(p± - p)j 



2)-pj 



j<0 
J>1 



(54) 



where p'^ w p— (tt— 2p)/Af andp^ w p+(7r+2p)/Af for a 



periodic chain of size M. The indices 



and 



cor- 



respond to phase slip and anti-phase slip, respectively, 
with the convention that the slip reduces the current. 
Schematically the saddle point and the metastable solu- 
tions are depicted in Fig. El Clearly the energy difference 




FIG. 6: Schematic representation of a metastable current car- 
rying state and an unstable saddle point solution. The arrows 
represent the superfluid phase at different sites of the lattice. 

between the two states is 

A£;_ =2JiV(2cosp -sinp(7r-2p)) (55) 

and 

AE+ = 2JN{2cosp +smp{TT + 2p)). (56) 

Correspondingly the decay rate to lower (higher) current 
state is proportional to 



Fzp CX e~'''^-^=f = g-2JAr/3(2cospT('rT2p)s 

In particular when p — > 7r/2 we have: 
F_ CX e-|^"''3(V2-p)='^ 



(57) 



(58) 
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In the one dimensional case it is also straightforward to 
evaluate the prefactor in the decay rate. We give details 



of such derivation in the Appendix |21 and quote only the 
final result here: 



J 



64 cos p jJN 
IT 



exp 



tt/2-p 



tanp - 



AN J 
T 



(cosp — (7r/2 — p) sinp) 



(59) 



64 cos p jJN 

1j 



exp 



7r/2 +p 



tanp 

r 



T 



(cosp + (7r/2 + p) smp) 



(60) 



Here r is a relaxation time of the condensate towards 
thermal equilibrium, which we will leave as a phenomeno- 
logical parameter. 

We can now compare the leading exponential terms 
for the two current decay mechanisms. Thus if we re- 
quire that the exponent in l|57|) is equal to the tunneling 
action computed in the previous section we can find when 
the two exponents coincide. It is convenient to introduce 
a characteristic temperature scale Tq, at which the en- 
ergy of the zero-point fluctuations is equal to the thermal 
energy of the corresponding classical system. Using Bo- 
goliubov's approximation one finds that 



T, 



2V2 



Q 



\fNJU. 



(61) 



We can now rewrite the expression for r_ at p close to 
Pc = 7r/2 in a more convenient form: 



64cosp / JN 



■ exp 



V27: UTjTq 3 



(62) 

Note that the exponent in the expression above coincides 
with that in (|5()|l when 



r*«0.21TQvW2^ 



(63) 



The temperature T* separates regimes of thermally acti- 
vated and quantum current decay. Thus if T < ther- 
mal phase slips are unimportant so that quantum tunnel- 
ing dominates the decay and vice versa. Note that unless 
p is very close to 7r/2, the crossover temperature T* is 
of the order of the characteristic Josephson energy Tq 
(or equivalently the sound velocity in the lattice units). 
Under present experimental conditions it is very easy to 
achieve T <C Tq and thus T < T* and therefore to ob- 
serve current damping due to quantum phase slips. 

In higher dimensions we can not find an explicit an- 
alytic expression for the energy of the metastable state. 
However, we can get an approximate result near the crit- 
ical current. Using again the idea that the transverse 
fluctuations can be treated in the continuum approxima- 
tion and expanding cos(p -I- (t>j+i — (pj) up to the third 



order in phases we can write the energy in the approxi- 
mate form: 



+ COs{p) -(t>jf 



pj+1 



,(64) 



where (x) is the nontrivial solution of the correspond- 
ing Euler-Lagrange equations vanishing at a; ^ oo. After 
rescaling (j) = cos(p) (f> and x = x^/2/ y/cosp wc find: 



i f U(^...-0.)^4(0..r-^.)^ 



.(65) 



Note that the integral in the expression above coincides 
with Sd~i up to a number . So using results H50|l . 
(|51|l . and (|58|l we immediately conclude that 



Si «1.3JAr(|_p)', 



~1Q JN[--p 



E3 ~ 35 JiV 



(f 



P 



(66) 
(67) 
(68) 



Correspondingly the exponents in the thermal and quan- 
tum decay rates become the same at a temperature 



T* « 0.44Tq\/V2 - 



P 



(69) 



both in two and three dimensions. This crossover tem- 
perature is about a factor of two higher than in the one 
dimensional case (I63II . If p is not too close to 7r/2, T* is 
again of the order of Tq and thus the quantum tunneling 
should be responsible for the current damping below the 
mean field transition at experimentally relevant temper- 
atures. 



13 



B. Current decay in the vicinity of the SF-MI 
phase transition. 

As we show in the Appendix IXI the quantum action in 
imaginary time takes the following form: 



S = C d'^zdx 



dV 

dx 



(70) 



Here z denotes the imaginary time and transverse coor- 
dinates which form a c?-dimensional space, a; is a longi- 
tudinal coordinate along the current. Note that in this 
section we measure coordinates in units of the coherence 
length. This is because we focus on the commensurate 
case and hence we are interested only in the superfluid 
regime r > (see Eq. In this case it is convenient to 

rescale x xj \pr to simplify the notations (in the orig- 
inal lattice units x is measured in the units of the corre- 
lation length ^) . The constant C depends on the original 
microscopic parameters of the underlying Hamiltonian. 
Within the variational ansatz described in Appendix ^ 
we findi2£: 



u 2(2d)''/2 



(1 



(71) 



where u = U/8dJN is the dimensionless interaction; 



^/2d(T^^' 



(72) 



In the case when thermal fluctuations are more impor- 
tant than the zero point motion, we are interested in the 
energy functional rather than the action: 



£^C' 



1 d'^-hdx 


dtp 


2 

-f 


dtp 


dz 




dx 



1 



\M' + 7;m\ (73) 



fluid transition, which is also of the Kosterlitz-Thouless 
class was studied in Ref. |51|. It was shown that the dis- 
sipation comes from unbinding of existing vortices and it 
does not have an exponential suppression. So while the 
mean field description near the quantum critical point in 
one dimension is questionable, we believe that it is justi- 
fied in two and especially in three spatial dimensions. 

Decay of superconducting current in the GL model was 
studied by several authors. In particular, the exponent 
characterizing the decay rate in the one dimensional case 
was studied in Ref. |30l | and t he p refactor setting the time 
scale was later found in Ref. [31| . In three dimensions at 
small currents the corresponding exponent was derived 
by Langer and Fisher—, where it was shown that the 
decay rate vanishes as exp(— C/p) as p — > 0. However, 
here we are interested in quite the opposite limit p pc, 
where the method used in that paper does not work. 

Let us start our analysis from the quantum decay. We 
first emphasize that quantum in this context means due 
to fluctuations beyond the saddle point approximation. 
Contrary to the Gross-Pitaveskii regime, where the clas- 
sical description is well controlled by the smallness of the 
ratio U/JN, there is no obvious small parameter here. 
The validity of the mean field description in this case 
can be checked a-posteriori by explicit computation of 
quantum corrections. The other comment we would like 
to make is that the parameters C and ^ entering the 
Ginzburg-Landau action are generally renormalized and 
deviate from the mean field results. 

To compute the tunneling action we need to find a 
bounce-solution of the classical equations of motion in 
imaginary time. Instead of using complex fields ipj we 
introduce two real fields rj and describing amplitude 
and phase fluctuations around the metastable minimum: 



where z comprises now d—1 transverse coordinates only. 
The value of the constant C" can be found within the 
mean field approximation (see Appendix IXI): 



C = JN- 



1 



'.{2d)<i/- 



2-d/2 



(74) 



Before proceeding we would like to point out that the 
mean field expressions for ^ and C are valid in the vicinity 
of the quantum phase transition at large spatial dimen- 
sions. For example, at d = 1 the superfiuid insulator 
transition belongs to the Kosterlitz-Thouless universal- 
ity class and it is characterized by proliferation of vor- 
tices^Si^. In particular, dissipation in two spatial dimen- 
sions in the vicinity of the thermal superfiuid-to-normal 



(75) 



Here we intentionally use another notation k for the con- 
densate momentum, because it is measured in the units 
of inverse correlation length ^. It is related to the usual 
momentum p in inverse lattice units as k = p^,. We can 
expect that close to the critical current both rj and (f> 
remain small and we can expand the action up to the cu- 
bic terms in these fields to find the correct asymptotical 
behavior of the instanton action at k ky. 



J 



J dzdx (dzri)'^ 



V3 



V^d^^+-r)\ (76) 
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It is easier to start the analysis of lf7S|l calculating the 
corresponding energy in d = 1 . This problem was already 
solved in Ref. for all values of fc. So we will use the 
asymptotic of their expression at k —^ kc to compare with 
our result and thus check the validity of our scheme. It 
is easy to see that upon the transformation: 



l-fc2 



(77) 



the cross term in 77 and (/) vanishes in the quadratic part 
of (|76() . Since the amplitude 77 mode remains gapped at 
k ~ kc it can be disregarded and the approximate energy 
(which is equivalent to the action (|76|l at d = 0) takes 
the form: 

dx^{dl^f+2V3{k,-k)(dxcbf~^{dx^)\ (78) 



Upon rescaling: 



3t 



3* 2\/kc — k 



4>-^\/kc-k 



(79) 



we obtain 



9 33(fc,-fc)5/2 / dx{dlcj>f + {dx<jyf~{dx<pf. (80) 



The last integral is just a number equal to the energy 
of the saddle point, which is a nontrivial solution of the 
Euler-Lagrange equations: 



-2^ + 2(/.'-30'2 = O, 



(81) 



where 



d(j)/dx. This equation has a simple solution 
1 



cosh'^ x/2 ' 
which after substitution into (|Sn|) gives: 



(82) 



(83) 



This result coincides with that derived in Ref. [3l| at A: ^ 
kc- Now we can proceed with a general d-dimensional 
case. As in the previous section the instanton action 
in d + 1 dimensions is equal to the barrier energy in d- 
dimcnsions. After performing the transformations H77() 
and ignoring the gapped amplitude mode 77 we find: 

d'^zdx\{dl,4>f + '^{d^<t>f + \{dl<t>f 

+ 2%/3(fcc - k){dxc^f - ^{d^c^f. (84) 
v3 

It is easy to see that the rescaling required to make the 
action independent of the momentum is: 



33 2Vfc~^ 



0-2" 



6(fc,-A:)- 
(85) 



Upon these transformations the first term in the action 
becomes irrelevant and in the leading order in fee — fc the 
action reads: 



Sd = 



2d 



-{k,~k)^ J dzdx{\/cl)f + {dl(f,f-{dxcl))\ 

.(86) 

where V — {dz,dx) is the gradient in d + 1 dimensions. 
For the classical energy, as we mentioned above, one has 
to substitute d ^ d — 1 in (|86|l . The important difference 
between the result for the continuum model 1)86(1 and the 
lattice result (|44|l is that the power of the fcc — fc in the 
prefactor in (|86|1 is smaller than that in H44() . Moreover 
in d = 3 the whole scaling breaks down suggesting the 
the instanton action becomes discontinuous at d = 3 near 
the Mott transition. 

We can evaluate the integral in ((86(1 using the vari- 
ational ansatz. For simplicity we will take a separable 
function 



{z,x) ^A{z) 



tanh ax 
cosh (3x ' 



(87) 



where a and (3 are the variational parameters and 
the function A{z) can be found solving remaining one- 
dimensional problem. With this choice it is easy to show 
that in one and two dimensions we obtain: 



aid ~ 0.32, pM ~ 0.53, ~ 73{k, - k)^ . 



a2d ~ 0.72, I3id ~ 1.08, S2d « 67(fcc - fc)^. 



(89) 



In three dimensions this variational ansatz gives ^3 — > 
0, which implies breaking of the scaling as fcc ^ fc. In- 
deed, the power of (fcc — fc) in ((86(1 becomes negative. 
However, it is unphysical to expect any divergence near 
the point of instability. This indicates that the orig- 
inal ansatz that the longitudinal coordinate scales as 
l/Vfcc — fc is not valid in this case and the tunneling ac- 
tion becomes momentum independent as fc — *■ fcc. To 
see that this is indeed the case and to estimate the ac- 
tual value of s in three dimensions we return to the the 
action ((84(1 without preforming rescaling 1(85(1 and apply 
the variational ansatz 1(87(1 . Then as fc ^ fcc we find: 



a3d ~ 0.53, Psd ~ 0.8, s^d ~ 125. 



(90) 



Using the mean field parameters C and ^ in l(70() and 
the results ((88(l - ((90l) we can rewrite the tunneling action 
in the following form: 



5.7 / ^ \3/2 



^2 
^3 



4.3.. 



(91) 

(92) 
(93) 



The equilibrium superfluid-insulator phase transition 
corresponds to ^ = 00. Notice that in one and two dimen- 
sions the tunneling action is always small as long as ^ ^ 1 
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and p is close to the critical momentum. This implies that 
at small currents the broadening of the nonequlibrium 
transition is very large. This is consistent with earlier 
numerical findingg^^ and recent experiments^^. In three 
dimensions, as we argued above, the tunneling action is 
discontinuous at the transition, therefore the mean field 
phase boundary is defined very well. So at d = 3 at zero 
temperature we can accurately define a relatively sharp 
phase transition between the current carrying superfluid 
and the insulator. 

Quite similar considerations apply to the current de- 
cay due to thermal fluctuations. The decay rate is deter- 
mined by the ratio of the activation energy and the 
temperature. In different spatial dimensions we get: 

E, ^ 1.3:1^ (l - VSp^f , (94) 

E, ^ 7.8^ (l - VSp^y , (95) 

Es- 8.6^ (l- VSp^y . (96) 

It is obvious that the thermal broadening is also strongest 
in one dimension. However, contrary to the quantum 
case, even at three dimensions close to the mean field 
transition the activation barrier vanishes continuously. 
Only in four dimensions and above we would be able to 
define a sharp phase boundary separating current carry- 
ing superfluid and insulating phases at finite tempera- 
ture. 

We would like to point out that the thermal decay 
is strongly suppressed at low temperatures T ^ JN/£^. 
Note that this condition is also necessary in order to ob- 
serve the equilibrium superfluid-insulator transition and 
thus can be satisfied experimentally. Another important 
point is that the action for the quantum phase slip tun- 
neling in one and two dimensions is never large near the 
mean field critical current. This implies that the quan- 
tum decay mechanism should be experimentally relevant 
at d = 1 and d — 2. This conclusion is similar to that 
we reached in the previous section when we discussed 
current decay at small interactions. 

It is also possible to make some qualitative statements 
beyond the mean field approximation in the vicinity of a 
quantum critical point separating equilibrium superfluid 
and insulating phases. Thus we can still expect that both 
the tunneling action and the thermal activation barrier 
vanish at 

P-^- (97) 

On the other hand, quite generally 1/^ ^ X'^, where j/ is a 
critical exponent'^^ and A is a tuning parameter across the 
transition, say deviation of the interaction U or hopping 
J from the critical point. In the one dimensional case 
v — oo [H^l (more precisely for the Kosterlits-Thouless 
transition ^ cx exp(6/VA), where b is some constant). 
In two and three dimensions the quantum critical point 



is characterized by the universality class of classical xy- 
model in one dimension higher and the corresponding 
critical exp onents are v w 0.67 a.t d — 2 and v — 0.5 at 
d — 3 "SS?!. We see that in three dimensions the mean 
field description gives the correct value of i'. Also, quite 
generally, we can argue that the action H7U|) and the en- 
ergy H73(l remain valid near the quantum critical point, 
but with constants C and C being renormalized: 

C cx (J, ^ ^d-4_ ^93) 

This in turn implies that near the quantum critical point 
A ^ 1 we get 

Ed ~ JNA'^Y'^ [B[iy - pf'^-"^, (99) 

where Ad, A'^, Bd, and B'^ are nonuniversal numbers. In 
the non-generic commensurate case, which we are mostly 
interested in here, z = 1. The scaling form above agrees 
with the mean field results {i/ = 1/2) obtained earlier. 
Despite quantitative difference between the correct and 
mean field scaling in one and two dimensions, the qual- 
itative features of the nonequilibrium phase transition 
discussed after and (|^ remain intact. 



V. DYNAMICS OF THE DECAY. INFLUENCE 
OF THE CONFINING POTENTIAL. 

A. Underdamped versus overdamped dynamics 

As we showed above, quantum or thermal fiuctuations 
lead to the broadening of the dynamical phase transi- 
tion. However, this does not imply that within a single 
experimental run a gradual current decay will be nec- 
essarily detected as the system is slowly tuned through 
the crossover region. The tunneling or the thermal ac- 
tivation times define a probability of generating a single 
phase slip. Once created the phase slip can either decay 
into phonon (Bogoliubov's) modes and bring the system 
to a next metastable minimum with a lower current, or 
this phase slip can trigger the current decay in the whole 
system. This situation is analogous to the motion of a 
particle on a tilted washboard potential with friction (see 
Fig.©. If the friction is large enough (or the tilt is small) 
then the particle, after it tunnels, will be stuck in the 
next minimum. On the other hand in the frictionless 
case a single tunneling event will cause accelerated mo- 
tion of the particle through the whole lattice. In a closed 
system these two regimes are well defined because the 
damping of the phase slip motion comes from the inter- 
nal degrees of freedom, which are completely described 
by the equations of motion. To see which of the regimes 
is realized in our systems within the classical thermal 
decay mechanism we numerically solve Gross-Pitaveskii 
equations of motion for an array of condensates. We start 
from a uniform current state. To have a current decay 
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FIG. 7: Possible motion of a particle after a tunneling event 
in a tilted washboard potential if the friction is small (left) 
and if the friction is large (right). Instead of changing the 
friction, one can vary the tilt. It is clear that reducing the tilt 
is similar to increasing the friction. 

we add small fluctuations to the initial values of the clas- 
sical fields. This is similar to starting from a thermal 
state. Since we cannot change the internal friction, in- 
stead we consider two different tilts. In Fig.|Slwe plot the 
computed current versus time for a one-dimensional ar- 
ray of M — 200 sites with periodic boundary conditions. 
The initial state is the noninteracting U = condensate 
with a given phase gradient p (specifically p — 2tt/5 and 
p — 7r/10) and unit hopping. It is clear from the figure 




Time (t) 



FIG. 8: Current (scaled to one at t = 0) versus time for a 
one dimensional periodic array of 200 sites with two different 
initial phase gradients. Here and in the following graphs time 
is dimensionless. Its units are set by the inverse units of cou- 
plings J and U. The evolution is determined solving equations 
of motion @ with constant hopping amplitude J — 1 and in- 
teraction increasing in time U = 0.01 tanh O.Olt for p = 2n/5 
and U = tanhO.Oli for p = tt/IO. To get the current decay 
we add small fluctuations to the initial values of the classical 
fields ■\l)j(t = 0). 

that we have an overdamped case for the smaller current 
case, where the phase slips occur one by one. On the 
other hand for the larger current a single phase slip gen- 
erates immediate current decay in the whole sample and 
this corresponds to the underdamped regime. We will 



not attempt here to find the precise boundary between 
the two scenarios, since it is not the purpose of our pa- 
per. We would like to emphasize that near the mean-field 
instability the system is always underdamped, while if 
the current decays at small p the motion is overdamped. 
These considerations agree with the experimental obser- 
vationsiiSi. We checked that the situation is similar in 
other spatial dimensions. Thus, if the current decays 
close to the mean field instability, then in a given ex- 
periment one will observe a sharp transition from the 
superfluid to the insulating regime. However, the precise 
point, where the current decays will depend on the de- 
tails of the experiment, for example on the rate of change 
of external parameters like tunneling and interaction or 
on the rate of change of the phase gradient p if the sys- 
tem is accelerated. On the other hand in the absence of 
any fluctuations the transition is very sharp and always 
occurs at p = tt/2. We can also perform a similar anal- 
ysis close to the SF-IN transition. The qualitative result 
that close to the modulational instability the phase-slip 
motion is underdamped remains correct. However, we 
should stress that in one and two dimensions broadening 
of the mean field transition is very strong and the actual 
decay may occur very far from the critical current. In 
this case an overdamped scenario should be realized. 

Unfortunately we cannot simulate the dynamics of the 
decay due to quantum tunneling. However, we would 
like to argue, that near the critical current the fate of 
quantum and thermal phase slips is identical. This is 
because the tunneling (activation) barrier is very narrow 
and the classically allowed motion after the tunneling 
event starts very close to the maximum of the barrier 
(see Fig.|7||. 

We have to make another important remark that if the 
motion of phase slips is underdamped then in a truly in- 
finite system the current state is always unstable. Indeed 
the probability of a phase slip is proportional to the sys- 
tem size M. If it causes the current decay in the whole 
system, then obviously a uniform current state cannot 
exist. However, in finite size systems these effects are 
not so crucial, because the decay probability depends ex- 
ponentially on the couplings and the current but only 
linearly in the system size. So the phase diagram plotted 
in Fig. O is quite robust to changes in M. 



B. Decay in a parabolic trap. 

If in addition to the optical lattice potential the con- 
densate is placed into a parabolic trap, then even at the 
classical (Gross-Pitaevskii) level the condensate momen- 
tum is not a conserved quantity. In a typical experi- 
ment, where the SF-IN transition is probed, the lattice 
potential is slowly ramped up resulting in hopping am- 
plitude going down. In Appendix |21 we show that in this 
case the amplitude of momentum oscillations increases 
in time: p oc (1/J(t))i/'* (see Eq. (EHJ). If we ig- 
nore completely the effects of quantum fiuctuations, then 



17 



for arbitrarily small initial displacement the condensate 
momentum will ultimately exceed the critical value of 
7r/2 and the condensate motion will become unstable. 
If this happens before the quantum fluctuations become 
significant, this effect will dominate the SF-IN transi- 
tion. In Fig. El we show the time evolution of the center 




FIG. 9: Center of mass momentum (in the units of tt per 
lattice site) versus time for a two-dimensional condensate in 
a parabolic trap with hopping amplitude decreasing in time: 
J{t) = 2.5 exp(— O.Olt). The other parameters are U = 1, 
number of atoms per central site A^'o — 1.5, strength of the 
confining potential k = 0.02. 

of mass momentum pcm = J2pP''^p/ J2p''^p of 
densate in a trap using Gutzwiller approximation. Here 
Up = ;(ajai) exp(ip (/ — j)) is the momentum distri- 
bution function. The condensate initially in equilibrium 
is given a small momentum boost. In agreement with 
our expectations the amplitude of momentum oscillations 
grows in time. We note that at the same time the conden- 
sate velocity v{t) cx J{t)p{t) decreases with time. This 
behavior continues until the momentum exceeds a critical 
value, where the condensate motion becomes chaotic. 

One can avoid complications related to the noncon- 
servation of the amplitude of momentum oscillations by 
tuning the interaction strength rather then the hopping 
amplitude. In this case one can directly probe the bound- 
ary separating stable and unstable motion of condensates 
at a given condensate momentum. Another important 
feature, which distinguishes trapped condensates from 
homogeneous systems is the spatial variation of the den- 
sity. Thus if the density profile is smooth enough, the 
condensate motion becomes unstable first near the edges 
where N k 1. In the center of the trap current decays 
at higher interactions. As we argued earlier in this sec- 
tion, in homogeneous systems the motion of phase slips 
is underdamped near the instability, i.e. a single phase 
slip triggers current decay in the whole system. We can 
expect the same to be true in a trap as long as the mean 
occupation number in the central site remains close 



to unity. On the other hand if A^o ^ 1 then it is in- 
tuitively clear that phase slips occurring near the edges 
can not destabilize the motion of the bulk of the con- 
densate, which is very far from the instability. To verify 
this reasoning numerically we again employ Gutzwiller 
approximation. In Fig. IIUI we plot the momentum oscil- 
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FIG. 10: Time dependence of the condensate momentum in 
a two-dimensional harmonic trap with different filling factors 
per central site. 

lations of a two-dimensional condensate versus time. We 
set the hopping amplitude J = 1/4 while increasing the 
interaction linearly in time: U{t) = O.Olt. The simula- 
tions are performed on a lattice of dimensions 120 x 60 
with global trapping potential V{jx,jy) = 0.01(j^ + Jy)- 
We consider two different filling factors in the central 
site A'o = 1.5 and A^o = 3. It is obvious that the onset 
of instability in both cases occurs at roughly the same 
interaction strength (which is close to the uniform re- 
sult at filling factor A^ = 1). However, while the motion 
becomes chaotic very fast for A^o = 1-5, there is a very 
gradual decay of momentum oscillations at larger filling 
Nq = 3. In agreement with our expectations this indi- 
cates that an overdamped current decay is realized in the 
latter case. 



VI. LOSS OF COHERENCE IN THE 
NONEQULIBRIUM PHASE TRANSITION. 

The superfiuid to Mott insulator transition at equi- 
librium, is a continuous quantum phase transition. As 
such, it is expected to be reversible. That is, if we tune 
through the transition and then back to the initial state 
slowly enough, we would recover arbitrarily large frac- 
tions of the superfiuid density^. 

With finite current the situation is markedly different. 
We envision the following experimental procedure. (1) 
The condensate is boosted to small but finite velocity 



18 



where 




FIG. 11: Possible experimental sequence to observe a super- 
fluid insulator transition in a moving condensate. 



(dipole oscillation in a trap). (2) An optical lattice is 
turned on adiabatically, until the system passes the line 
of instability (See Fig. [TTll . and then slowly turned off. 
Finally the atoms are released from the trap, and their 
final momentum distribution is measured and compared 
to the initial one. 

The current carrying state has a finite energy, which is 
released into the system when the current decays. Since 
the system is non integrable, it is reasonable to assume 
that the energy associated with the supercurrent, is re- 
leased in the form of incoherent excitations, which even- 
tually thermalize. If now the system parameters are 
changed sufficiently slowly, the subsequent evolution of 
the system is adiabatic conserving the entropy of the sys- 
tem. Thus, the entropy change following the current de- 
cay, may be used to obtain the depletion of the superfluid 
in the final state of the system. We would like to empha- 
size that these assumptions apply if the parameters of the 
system change not too fast, so that the current decays in 
a quasi static regime. Otherwise the current decay and 
the entropy release are not governed by the properties 
of the transition point and the actual time dependent 
problem has to be solved. 

The temperature of the equilibrium state reached fol- 
lowing the current decay, and before the system param- 
eters had a chance to change appreciably may be calcu- 
lated by equating the energy of the superflow, prior to 
its decay, with the internal energy of the system in the 
new thermal equilibrium: 



Se{p) 



-E- 



(100) 



The low energy excitations in the superfluid state are lin- 
early dispersing sound modes with Wk ~ cfc. Substituting 
this in (|100|l . we can solve for the temperature to obtain: 



(101) 



A,= 



nadlQil + d) 



(102) 



Vld is the surface area of the d-dimensional unit sphere, 
and C stands for the Riemann zeta function. Accordingly, 
the entropy of this thermal state is given by: 



S 



Hp) 

T 



Se 



1 



(103) 

As argued above, we can use this entropy to calculate the 
condensate depletion in the final state. 

Let us apply this procedure assuming that the current 
decays via the instability in the vicinity of the Mott tran- 
sition. Though we do the calculation for all dimensions, 
one should note that such a scenario is particularly rel- 
evant for a three dimensional optical lattice. According 
to section in lower dimensions the current will most 
probably have decayed due to fluctuations before reach- 
ing the instability. 

The energy per site of a state near the Mott transi- 
tion according to the Ginsburg-Landau model, (see Ap- 
pendix is given by: 



JN 
2duV 



d'^2 



(104) 



Substituting the field corresponding to the current car- 



rying state: -i/i = ^ — 
5e = e{p) ~ e(0) 



we obtain 



(,-2 1 2. 2 



(105) 



Recall that our proposed experiment maintains constant 
p and changes the dimensionless interaction, and hence 
also ^, by increasing the lattice intensity. The current 
decays when the instability is reached, i.e. when = 
^jT^ — 3p^, at which point the energy per site is: 



6e 



JN5 4 
'2d 2^ 



(106) 



Thus the energy released following decay via the insta- 
bility at phase gradient p is oc p^. Using (|103|) and the 
sound velocity near the transition c — 2JN^/2d we get 
the increase of entropy per site: 



S 



1 ( 



d-Ai V4(2c/)3/2 



(107) 



In three dimensions, this gives S-^d ~ O.lGp'^. With 
such a small increase of entropy we anticipate that the 
irreversibility, as manifest in the unrecovered condensate 
fraction, would also be small, for low initial currents. 
Perhaps the simplest way to see this is to slowly reduce 
the lattice intensity until the atoms are very weakly in- 
teracting before releasing to measure the momentum dis- 
tribution. In this case the elementary excitations have a 



19 



quadratic dispersion uj^ = ak'^ . In general this assump- 
tion is not necessary and one can use the full Bogoliubov 
spectrum. However, qualitatively the result remains the 
same. A nice feature of the quadratic dispersion, is that 
the thermal depletion is simply related to the entropy, 
which is given by: 



S = 



Where 



d + 2\ E 



d J T (27r)" 



d + 2 



ji\ d/2 

Id (108) 

a J 



Id = 



dx- 



(109) 



The thermal depiction on the other hand is 
q 



t r dxx"-^^ 



(110) 



where xq depends on the small momentum cutoff deter- 
mined by the system size (i.e. xq ^ L^Ja/T). In three 
dimensions the last integral is convergent and we get: 



,3d 



2C(3/2) 
5C(5/2) 



S w 0.785' 



(111) 



Thus, in three dimensions the number of excited quasi- 
particles, or condensate depletion, is a direct measure of 
the entropy 

In one and two dimensions on the other hand, the inte- 
gral has an infrared divergence. Of course this is simply a 
restatement of the well known fact that a true condensate 
in free space at finite temperature does not exist below 
three dimensions. In practice, the divergence is cut off by 
the system size and one can define the notion of a qua- 
sicondensate. In two dimensions where the divergence is 
only logarithmic: 
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And finally in one dimension 



Id 



(112) 



(113) 



Instead of the system size the cutoff may come from the 
finite momentum resolution (Ap) of the experimental ap- 
paratus. In general the momentum cutoff will be deter- 
mined by the minimum of L and 1 / Ap. We see that below 
three dimensions condensate depletion due to thermal de- 
cay at small currents is more pronounced and easier to 
detect than at d = 3. 

The situation is different if the current decays at 
smaller interaction strengths (smaller lattice intensity) 
at large currents. The energy of the current state may 
then be calculated from the Gross-Pitaevskii energy func- 
tional: 



(ij) » 



(114) 



Substituting VNe'P''^ for ipj we find Se = 4 JiV sin^(p/2). 
Together with Eq. ((TUSjl . and using the sound velocity 
mtrop; 

(—\ sin|, (115) 



c — \/2UJN we get the entropy increase per site 



S,,^2.2y—J sm' 2, 



S,d ^ 2.2 y^— J sm" ^. 



(117) 



Evidently the irreversibility is stronger and easier to de- 
tect if the current decays at smaller interaction strengths. 



VII. EXACT RESULTS IN SMALL SYSTEMS. 

In this section we present results of exact dynamics in 
small one-dimensional systems. We will assume that we 
have a periodic one-dimensional array of M sites contain- 
ing N particles. At i = we assume that interactions are 
absent and the system is placed in the uniform current 
state, described by the wave function: 



1 



N 



|0), 



(118) 



where |0) is the no-particle vacuum. This state is an 
exact eigenstate of the noninteracting system. Now we 
slowly turn on interactions driving the system closer to 
the insulating regime and then gradually turn them off. 
The latter step is necessary to check whether we have re- 
versible dynamics or irreversible current decay. Although 
in real experiments it is rather tunneling not the inter- 
action which is changed in time, this does not make any 
qualitative difference in uniform systems. If the interac- 
tion is ramped up infinitesimally slowly, then any finite 
system will remain within a particular energy eigen-state 
and the evolution will be always reversible. However, the 
energy splitting between the many body levels decreases 
exponentially with the system size and the number of 
particles. So practically even in relatively small systems, 
one can study dynamics, which is slow with respect to the 
characteristic time scales (like period of Josephson oscil- 
lations), but very fast with respect to the inverse many 
body energy spacing. To make this point more transpar- 
ent we plot in Fig.^lthe energy spectrum versus U / J for 
a particular system of a periodic array consisting of six 
sites and containing six particles. The size of the Hilbert 
space here is already quite big: TV = ll!/(6!5!) = 462. 
It is obvious from the figure that while the ground state 
is well separated from the continuum at all interaction 
strengths, the excited states experience many level cross- 
ings. So it is nearly impossible to be in or close to the 
adiabatic regime unless the system is in the ground state. 

Let us assume that the Hamitonian is described by , 
where the interaction strength changes in time according 
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FIG. 12: First one hundred energy levels of the spectrum 
in a periodic array of six sites containing six particles (top) 
and zoom of the energy spectrum around the current car- 
rying state with p = n/Q (bottom). It is clear that while 
the ground state is well separated from the continuum for all 
interaction strengths, the excited many body states are com- 
pletely mixed. So it is nearly impossible to be in the adiabatic 
regime unless the system is initially in the ground state. 
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FIG. 13: Dependence of current on time for the system of 
eight sites with two particles per site. The hopping amplitude 
is equal to unity and is independent of time. The interaction 
changes in time according to Eq. 111911 with C/o=2, 5 = 0.02, 
and T = 200. 
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U{t) = Uottmh5t ia.nh5{T ~t), (119) 

where 5 is the adiabaticity parameter, T is the duration 
of the time evolution and is a prefactor setting the 
energy scale. We also assume that the hopping J is equal 
to unity and does not change in time. 

In Fig. ll3l we plot the current versus time for the array 
of eight sites with two particles per site. The two curves 
correspond to initial phase gradient of 7r/4 and 7r/2 per 
bond. The smaller current decays when the interaction 
becomes sufficiently large, in agreement with that this 
state is metastable, while the 7r/2 state decays almost 
instantly since it is classically unstable. In both cases 
the decay is clearly completely irreversible. To make the 
final point more transparent we also plot momentum dis- 
tribution as a function of time for the same parameters in 



FIG. 14: Occupation of momentum states as a function of 
time for different current carrying states. Solid lines corre- 
spond to the occupation of momentum p equal to the initial 
phase gradient in the system. The dashed lines are the occu- 
pations of the zero-momentum state. The parameters are the 
same as in Fig. 1131 



Fig. 1141 The curve corresponding to a zero current state 
shows reversible behavior, suggesting that the interac- 
tion indeed changes slowly enough so that the system 
is in the adiabatic limit. Note, that even at the peak 
interaction strength, the system is still far from the insu- 
lating state, as evident from the very small depletion of 
the p = state. Nevertheless, because the smallest phase 
difference per site achievable in an array of size eight is 
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still quite large {p = 7r/4), this interaction is sufficient to 
drive current decay. Indeed the curve corresponding to 
initial phase gradient of 7r/4 clearly displays metastable 
behavior, with current decaying after some delay. In the 
new steady state reached, occupation of zero momentum 
builds up suggesting that it corresponds to a thermal 
distribution with no macroscopic current, but with some 
residual phase coherence. The state with phase gradient 
p — tt/2, by contrast, seems to decay into a high temper- 
ature state, without visible phase correlations. Indeed, 
occupation of momenta zero and tt/2 are equal to unity 
as expected if the phases were completely random. It is 
peculiar that there are only very weak fluctuations of the 
momentum distribution in this state. This should be con- 
trasted with the rather large fluctuations seen following 
decay of the ir/A current state. 

Finally Condensates sustaining phase gradients p > 
7r/2 are classically unstable, and are expected to decay 
rapidly, resulting in even higher temperature than for the 
7r/2 state. The p = tt state is an interesting exception to 
this rule. As argued in>^ji^, this state evolves into a max- 
imally entangled Schrodinger cat-like state, which is ro- 
bust to weak perturbations in small systems. Physically 
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FIG. 15: Occupations of the k = and k = n states for the 
state with initial phase gradient p = n. Here the interaction 
changes in time according to ' l|119|l with Uo = l, S = 0.02, and 
T = 200. 

this happens because the tt state in the noninteracting 
case is the most excited state in the system. In the ab- 
sence of energy relaxation the system remains in the most 
excited state under slow perturbations in the same way 
as it does in the ground state. And therefore we can ex- 
pect reversible behavior of the phase coherence. We plot 
the corresponding momentum distribution for this state 
in Fig. El Although for the intermediate times the evo- 
lution looks completely chaotic, once the interaction is 
reduced back to zero the momentum occupancy a,t k = tt 



reaches almost the maximum possible initial value sug- 
gesting only a small amount of excitations in the system. 

Unfortunately doing exact numerical simulations we 
are quite limited by the total system sizes and the number 
of particles. Also we can address only one-dimensional 
systems. Therefore we can not directly test the phase 
diagram and decay rates derived in the previous sections. 
Nevertheless, we would like to point that the numerical 
results are in excellent qualitative agreement with our 
predictions. 



VIII. SUMMARY AND EXPERIMENTAL 
IMPLICATIONS. 

In summary we emphasize two important predictions 
of this work. In Sec. lIIII we derived a mean field phase dia- 
gram for the stability of a moving condensate. In Sec. II VI 
we investigated the effect of quantum and thermal fluc- 
tuations, which lead to broadening of the nonequilibrium 
transition. 

We showed that the mean field transition interpolates 
between the classical modulational instability, which oc- 
curs at phase gradient p = 7r/2 deep in the superfluid 
regime and the equilibrium {p = 0) transition to the 
Mott insulating phase at strong interactions. The dy- 
namical transition is of first order (i.e. irreversible) at 
any nonzero current, contrary to the second order tran- 
sition at equilibrium. Thus, if one starts from a current 
state and slowly drives the system towards the insulating 
regime, e.g. by ramping up the lattice potential, then af- 
ter crossing the transition boundary the current decays 
irreversibly, releasing the energy of the coherent motion 
into heat. Plotting the location of the nonequilibrium 
transition as a function of the current and extrapolating 
the curve into the static regime p = is a way to accu- 
rately determine the position of the equilibrium SF-IN 
transition. 

The mean field theory does not take into account quan- 
tum tunneling and thermal activation of phase slips. 
These induce decay of supercurrent, even before the clas- 
sical equations of motion become unstable. We calcu- 
lated the asymptotic decay rates near the mean field in- 
stability in two regimes: (i) deep in the superfluid regime 
and (ii) close to the equilibrium Mott transition. 

In a three dimensional optical lattice the broadening 
of the transition due to these effects is found to be small 
in all cases. In particular we find a discontinuity in the 
current decay rate across the mean field transition for 
small currents at zero temperature (close to the equilib- 
rium Mott transition). Thus the dynamical transition 
survives the effect of quantum fluctuations in this case. 
We predict that a sharp dynamical superfluid-insulator 
transition would be seen at small currents at a critical 
interaction strength (or lattice depth). 

In one and two dimensions, on the other hand, quan- 
tum and thermal phase slips lead to substantial broad- 
ening of the transition, especially when the average site 
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occupation '--^ 1. Then we expect the current to decay 
weU before the dynamical instabihty is reached. Indeed, 
in a recent experiment^ strong damping was detected 
at currents much smaller than that given by the Gross- 
Pitaevskii modulational instability. In addition the ob- 
served dependence of the damping rate on the lattice 
depth potential was very smooth. This is in line with 
our prediction of a the large broadening of the mean-field 
transition by quantum fluctuations and should be con- 
trasted with the Gross-Pitaveskii predictions of a sharp 
transition. It is also consistent, with earlier numerical 
results by one of usM'. 

The experimental results^, do not by themselves prove 
that quantum rather than thermal fluctuations are re- 
sponsible for the observed damping. Here we estimate 
the crossover temperature from thermal to quantum 
dominated decay to be of order the Josephson frequency 
T* w V JU N /ks- The experiment is therefore likely to 
be dominated by quantum phase slips. To verify this 
conclusions, one could measure e.g. the damping as a 
function of temperature and observe a saturation of the 
rate around T* . 

One perhaps surprising experimental observation was 
that in the overdamped regime (i.e. at high optical lat- 
tice potential) the condensate was essentially localized 
in a tilted lattice, while it still exhibited sharp coher- 
ence peaks^. A possible explanation is the effect of 
the inhomogeneous density in a harmonic trap, which we 
discussed in Sec. Indeed, in the overdamped regime 
one expects no suppression of phase slip tunneling at the 
edges of the condensate, where N 1. This implies that 
the phase at the edges fluctuates very wildly and the edge 
atoms are localized. On the other hand, in the middle of 
the condensate, where the mean number of particles per 
site is larger, the system is far from the mean field transi- 
tion and phase slips are relatively costly. As a result, the 
edges of the condensate create an effective potential bar- 
rier stopping the motion of the rest of the system, which 
retains phase coherence. 

Our results are consistent with another recent exper- 
iment, where the superfluid decay was measured as a 
function of the condensate velocity in a one-dimensional 
optical lattice^''. There the average number of bosons per 
site was large and quantum fluctuations negligible. At 
low temperature a dynamical localization transition con- 
sistent with Gross-Pitaveskii predictions was observed. 
However, at relatively high temperatures, the motion be- 
came unstable at much lower quasi-momentum. This 
observation qualitatively agrees with the decay mecha- 
nism due to thermal phase-slips considered in the present 
work. 

We also presented exact numerical simulation of small 
one-dimensional systems. These were in qualitative 
agreement with the physical picture discussed in the pa- 
per. For example we demonstrated that at nonzero cur- 
rent the transition is irreversible. We also find that in 
a periodic chain of eight sites, the current state with 
p = 7r/4 decays only at some finite interaction strength. 



while at p = 7r/2 the decay occurs almost instantly. An 
important exception is the case with p = tt, where the 
evolution is reversible (see also Ref. |53). A quantita- 
tive comparison of the exact numerical results with our 
predictions is not possible because of strong finite size 
effects in the exact simulations. 

Note added - After this paper has been submitted two 
preprints appeared^SiS, which address the experiment by 
Fertig et aU^. There, the damping was attributed to sin- 
gle particle Bloch oscillations in the free fermion repre- 
sentation of the bosons in the limit of strong interactions. 
This effect can also be understood in the Boson language. 
If the number of particles in the trap is small, the system 
reaches the impenetrable boson regime while not yet in- 
sulating. Indeed one can easily transfer a boson from a 
filled to an empty site near the edge of the system. If the 
tunneling amplitude (J) is larger than the single particle 
energy near the edge of the cloud J > kN^ /8 then the 
created hole is delocalized through the whole system and 
the state is not insulating. Here kj'^/2 is the confining 
potential of the trap and Nt is the total number of par- 
ticles, which equals the system size in the fermionized 
regime. The hard core constraint in turn requires that 
U ^ J. Therefore if C/ ^ kN^ and we gradually de- 
crease J then indeed the system first goes to fermionized 
delocalized regime U ^ J ^ kN"^ and only if J is de- 
creased further it becomes localized. On the other hand 
if the total number of bosons is large U <C kN"^ , which is 
the case close to the thermodynamic limit, then the edge 
excitations in the fermionized regime are always localized 
and thus unimportant for the macroscopic properties of 
the system. 

We emphasize that if the first (i.e. small particle num- 
ber scenario) is realized, then after the trap minimum is 
displaced, the particles will essentially undergo Bloch os- 
cillations with different frequencies resulting in a damped 
center of mass motion^-. In this scenario there is no 
real energy relaxation of the center of mass and it will 
saturate at a displaced position^. In the second case 
U <^ kNf the system undergoes a Mott transition when 
the impenetrable regime is reached and the edge bosonic 
excitations can result only in a tiny center of mass dis- 
placement vanishing in the thermodynamic limit. The 
damping prior to the Mott transition in this case occurs 
via the mechanisms discussed in this paper. I.e., the cur- 
rent decay is irreversible and results in energy relaxation 
of the center of mass motion, which will eventually slide 
to the minimum of the trap. This seems to be the case 
realized in the experiment of Ref. ^2/i\ . 

We also mention that the single particle Bloch physics 
will dominate the decay mechanisms studied here if the 
Bloch oscillation, which frequency is equal to the single 
particle energy separation between the nearest sites due 
to external potential, is longer than the Josephson os- 
cillation. This condition is satisfied by the experimental 
systems of interes t "'"'^inl . 
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APPENDIX A: DERIVATION OF A 
GINZBURG-LANDAU ACTION NEAR THE 
SUPRFLUID INSULATOR TRANSITION 

A derivation of a Ginzburg-Landau action near the su- 
perfluid insulator transition was outlined in Ref. |39| |. 



J 



For convenience we present the full derivation in this ap- 
pendix. 

Let us choose the energy of a single site with integer 
N atoms, as the zero of energy. Then the Hamiltonian 
of the boson Hubbard model (QJ assumes the form: 



H 



U 



H.c.) + ^ - (n, - Nf ~ ii{n, - N) 

(Al) 

Close to the superfluid insulator transition, the particle 
number fluctuation is small. It is then possible to con- 
sider a subspace allowing only occupations of — 1, iV 
and A^ -t- 1 atoms per site. This reduced Hilbert space 
is conveniently described by a (over-complete) basis of 
product states: 



\n) [cos(0,/2) I A^) ^. + e"'^ sin(0,/2) (e'^^ sin(x/2) | A + 1 ) ^. + g-^'^^' cos(x/2) | A^ - 1 ) ^.)] . (A2) 



We shall use these states to construct a path integral for the evolution operator. The derivation can be carried out 
for arbitrary A, but for simplicity of presentation we take A^ >> 1. The first step is to prove a resolution of identity: 



pTT I'll plTT /'7r/2 

dd dx dip dr]M{n)\n) {n \ =i, 

Jo Jo Jo J --it/2 



(A3) 



We shall now find a suitable integration measure M{9), which is a function of 9 only. Substituting AA{6) in (|A3|) . we 
can integrate over r], ip and Xi which kills off the cross terms so that (|A3|) reduces to: 



I=2tt^ [ ddM{e)\ cos^-\N) {N\ +-'. 
Jo [ 2 2 



sin^ - lA^ + l) (A + l| + lA^-l) (A-1 



(A4) 



r 



The measure Ai (9) must enforce the identity between the 
diagonal matrix elements, so that: 



/" 

^0 



d9Mi9) I cos^ ^ - ^sin^- ) =0 



or equivalently 



(i6'A^(0)(l -I- 3 cos 61) 



(A5) 



(A6) 

where y = cos 9. This requirement is satisfied by A4{9) = 
C cos 9(3 cos 9 — 1), since it ensures that the integrand 
is an antisymmetric function of cos 9. The constant 
C = TT^'^ is determined from normalization. Since we 
are interested in the vicinity of the transition at f? = 0, 



where the measure A4 changes slowly, it is safe to replace 
it with a constant. 

It is also straightforward to calculate the Berry phase 

(1] I ^ I r!) = j^sin2(0,/2)(7?, - cosx,^,) = -«T(t). 

(A7) 

Equations ljA3p and ljA7|l are the necessary ingredients 
for the path integral representation of the evolution op- 
erator: 

U{t)^ Jvncyipi^iJ^dt'[T{t')-n{t')]y (A8) 
where the classical Hamiltonian is given by 
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,U 



TL ^ {fl\H\fl) = ^ sin^(6'i/2)(— - /icosxi) - 2JN'^piPj acj cos{rij - r]^ + tpi - tpj) + SiSj cos{rij - rjt + tpj - ipi) 



+CtSj cos(?7i + rij + (pj - ipi) + s^Cj cos{r], + tj^ + tpi - tpj) 



(A9) 



Here we introduced the notations q = cos(xi/2), = sin(xj/2) and pi = isin^^. It is important to note, that the 
dynamics defined by (jA7|) and ljA8|l . consists of two pairs of conjugate variables. The average offset from integer 
density is conjugate to the phase p^ while the second moment (i.e. the number fluctuation) is conjugate to 77. 

At fillings close to integer, the minimum classical energy is reached in a uniform state with x close to 7r/2 and 
77 = 0. We therefore expand the action to leading order in tJi = 7r/2 — Xi and rn- In addition we expand up to quartic 
order in p, and anticipating a diverging length scale at the transition, take the continuum limit of the action via a 
gradient expansion 



dt'd'^x [p^a{ip + m) - P^V ~ 2J7V {{Vpf + p^{Vipf) + AJNd{l - u)p^ - AJNd ■ up'^ - 2JNd ■ p^{2rf + ctVS)] , 

(AlO) 

where u — U/8JNd. We can now integrate over the gaussian fields 77 and a to obtain a Ginzburg-Landau (GL) action: 



S 



dt'd'^x 
1 



1 



AJNd 



AJNd 



dt'd 



(p2 + p^{p + fi)) ~ 2JN {{Vpf + p^{Vpf) - 4JNd{l - u)p^ + AJNd ■ up^ 



\{dt + ip)i^\'^ - (2JiV)22d|VV'P + {UNdf{l - u)\il;\'^ ~ {AJNdfu\tl)\'' 



(All) 



where tp = pe"^. Note, that in Eq. I|A10|I we left out 
terms of the form p'^{Vr])'^ and p^(Vct)^. After integrat- 
ing over 77 and ct, these would lead to irrelevant high order 
derivatives in the GL action. We can identify a sound ve- 
locity c = 2J7V\/2d, where the lattice constant is set to 
1 = 1. If we now make the transformation t ct (i.e 
measure time in units of l/c), and ?/' 'tpVTdu, the GL 
action assumes the form: 

S= - [ dtd''x\{dt+ip)^\^-\V^\'+r\i;\^-l\ij\\ 
a J 2 

(A12) 

where r = 2d{l — u) and a = 2u{2d)^'^. In the super- 
fluid phase r = <^~^, where f is the mean field coherence 
length. Note that the action is Lorcntz invariant only at 
commensurate filling, where by our choice of the zero of 
energy for (|Aip . ^ = 0. This is due to the particle hole 
symmetry, which is present only at commensurate filling. 

From the action ljA12p , we can find the deviation from 
commensurate filling, given by the conserved charge: 



5S 



Sfj,(t) ai 

(A13) 

In ijHJ, we chose the zero of energy such that the chemi- 
cal potential is /i = at commensurate filling N . Indeed 
if we substitute /i = in (jA13|l and a time indepen- 
dent order parameter we get Q = as required. How- 
ever, we note that the time derivative always appears in a 
gauge invariant combination with the chemical potential. 
Therefore the action l|A12|l and the "charge" l|A13|l . dot 
not depend on the particular choice of the zero of energy. 



In particular for calculating the dynamics it would prove 
convenient to use a different gauge, applying the transfor- 
mation ?/; — > T/je"^*-'-*, /i — > /i — (/) with = pt. This elimi- 
nates p from the action, at the expense of imposing on the 
order parameter an aditional time dependent phase. The 
two gauges coincide at the commensurate point where 
/i = and there is no time dependent phase. At in- 
commensurate filling, though the action seems Lorentz 
invariant in the new gauge, the physics is clearly not, 
due to the imposed time dependent phase 

We also note that in any gauge we can trace back 
the density parameter a = 7r/2 — x, appearing in the 
Gutzwiller states (|A2ll . In mean field theory, integrat- 
ing out a in (jA10|) . simply enforces the identity a = 
y^8/d{ip + p). A small incommensurate filling is then 
given by ap^. 

The Euler-Lagrange equations derived from the action 
in the new gauge (where p is eliminated from the ac- 
tion), are given by Eq. Hll|l . which we reproduce here for 
completeness 



(A14) 



In this gauge the density offset from integer filling is set 
exclusively by the initial conditions for ip, and and 
given by 



Q = 



1 



d'^x{'il;*ip — ipip*). 



(A15) 



It is easily verified that (|A15(I , is a conserved quantity in 
the equations of motion ljA14|l . This gauge choice, is thus 
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analogous to the canonical ensemble, where the particle 
number is fixed and is automatically conserved by the 
dynamics at all later times. 

Before concluding this appendix let us make several 
notes. First, to obtain the action l(7n|) from Eq. HA12|I . 
one has to rotate to imaginary time t — ixo, rescale 
length X ^ (i.e measure length in units of the coher- 
ence length), and the order parameter tp £,~^ip. Then 
in the "canonical gauge" the action assumes the form: 



2u(2d)3/2 



(A16) 



which is identical to Eq. (|7()|l . Second, to obtain the 
classical energy from (jA12ll in the original static gauge, 
we simply set time independent fields and multiply by 
the energy unit c/l: 



,(A17) 

The average density offset form integer filling is then 
given by: 



an a 



(A18) 



which agrees with with Eq. (|A13p evaluated in the static 
gauge, as it should. 

Another note we should make is that for low filling 
factors TV ~ 1, the general form of the action remains 
the same, in the vicinity of the Mott transition. However 
in that case, the expressions for ^, c, and the prefactor 
of the action are more complicated. 



Pitaveskii equations: 



^ I V.IV, (B2) 



Here Cj(t) is the Langevin noise term; ujj and oj^ are 
the eigenvalues of the excitations of Eq. H53|) around the 
saddle-point and metastable states respectively. Note 
that both states have one zero eigenvalue due to global 
phase U{\) symmetry and we ignore them. The sad- 
dle point spectrum also has one imaginary eigenvalue 
(wo) corresponding to an unstable solution of the lin- 
earized equations of motion. An extra prefactor w^^j 
comes because in the saddle-point configuration the num- 
ber of solutions with real u is smaller by one than in the 
metastable state. For simplicity we assume L to be even. 
Because of the absence of continuous translational sym- 
metry, there is no second zero eigenvalue for the saddle 
point (compare with Ref. [sJl). The eigenfunctions of 
small fluctuations around the metastable state are plain- 
waves. The corresponding spectrum thus reads: 



(B3) 



.L - 1 



u)n = 2y^2cosp sinfc„/2, 

where fc„ — 21:71/ L is the momentum, n — 0, 1, , 
is an integer and L is the size of the chain. 

The saddle point solution with a single phase slip has 
scattering states and a bound state. The latter is de- 
scribed by the eigen-function: 



A 



J> 1 
J<0. 



(B4) 



APPENDIX B: DERIVATION OF THE 
PREFACTOR OF THE CURRENT DECAY IN A 
ONE-DIMENSIONAL LATTICE 



Following a general theory developed in Ref. [45j , in the 
one dimensional case, we can find the transition rates per 
site r_ and r+. These can be written as: 




where r is the phenomcnological coupling to the bath 
degrees of freedom, defined in terms of dissipative Gross- 



Substituting ljB4p into the linearized version of H53() we 
find K — hii and ujq = Ai •\/2/3 cosp, where we ignored a 
small discrepancy between p and p' . 

To find the scattering states we shall seek solutions in 
the form: 

5(t)j = Ae'''^ + Be-'^i (B5) 

for J = 1, 2, ... L — 1 and 5(j)Q = d(j)is[, where S(j)j is the 
small deviation from the saddle point solution (|54|l . The 
system of secular equations determining the wave vectors 
kn reads: 



^ (1 -I- e^''^ - 2e''') +B{1 + c-'''^ - 2c-'''') = (B6) 
Ae''' (1 + e'''^ - 2e*''(^"^') + Be-''' (l + e-'''^ ~ 2e-'''^^-^'>^ ^ 0. (B7) 

I 

A nontrivial solution of the above system exists for k satisfying the following equation: 

kL k 
tan — = 2 tan - (B8) 
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Introducing the phase-shift kn = 2Trn/L + 5n/ L we find 

In the Umit L — > oo we get the approximate solutions for 
the scattering phaseshifts: 



2 arctan ( 2 tan 



The energies of the scattering states are equal to 



i>'n — 2\/2 cosp' sin 



L 



5^ 
2L 



(BIO) 



(Bll) 



Now it is straightforward to find the ratio of the products 
of all eigenvalues at the saddle point and the metastable 
state in the limit i — > oo: 



n 



3e 



± tan ; 



(B12) 



Substituting this into the general expression HB1() we de- 
rive (P|l and P?)). 



APPENDIX C: GROSS-PITAEVSKII DYNAMICS 

OF A LATTICE CONDENSATE IN A 
PARABOLIC TRAP UNDER SLOW CHANGE OF 
THE HOPPING AMPLITUDE. 

In the superfluid regime the most significant effect of 
the optical trap on the motion of the system is the mod- 
ification of the effective mass of the particles. So, let us 
analyze a condensate with a time dependent mass, mov- 
ing in a parabolic trap. Since we are considering a purely 
classical effect we can use Gross-Pitaveskii equations (for 
simplicity we restrict the analysis to one dimension: 



9V 



2m OX'' 2 



(CI) 



In the optical lattice the underlying equations are 



- J(Vj+i + ^j-i) + ttjV + C/|^2|V-, (C2) 



In the weak tunneling regime the equations HC1|I and 
(|C2|I are equivalent provided that m = 1/(2 J), the lat- 
tice constant is equal to unity, and the phase gradient is 
small. If the interaction strength is not too small, i.e. the 
condensate is in the quantum-rotor limit UN ^ J , then 
the effect of quantum pressure is negligible and instead of 
(|C1|) one can use hydrodynamic equations of motion^. 
The bosonic field is then represented as: 



(C3) 



Keeping only the lowest orders of spatial derivatives of 
the density p instead of \CA\ we obtain: 



dp 
di 
dp 
dt 



-kx — 

)_d_ 

m dx 



1 dp dp 



rrJ' dx 
(PP)- 



U 



dx^ 



(C4) 
(C5) 



The stationary solution of equations ljC4p and (|C5|I yields 
an inverted parabola profile of the density: 



Poix) 



2" kx 



u 



(C6) 



where fi is the chemical potential. Let us now assume 
that the condensate undergoes small center of mass os- 
cillations. They can be excited by a small displacement 
of the trap minimum. Then it is easy to check that one 
can seek a solution in the form: 



p{x,t) 



k 2 
2 



f{t)x 



A 



dpjt) 
dt ' 



(C7) 



with the initial conditions: /(O) — kxQ, p{0) — 0, 
/i(0) — po ~ ^fcxQ , where xq is the initial displacement. 
Substituting l|C7|l into l)C4|l and l)C5|l one finds: 



f{t) = xok cos \ —t, p{t) = 
m 

_ m \2 

p{x,t)^— ^ ' - ' 



-Xovfcmsin \ — i, 

TO 



U 



(C8) 



Note that the interaction U never enters (|C8(I except for 
a trivial prefactor. This is so because the excited mode 
is related to the motion of the center of mass, while 
the shape of the condensate cloud does not change in 
time. Moreover the existence of such a center of mass or 
Galilean mode is the property of the equation ljCl|) itself. 
Indeed, it is easy to check that if ^"0(2;,^) is the solution 
of (Ell) then 

7/;(a;, t) = - so(t), i)e'P(*^^e-5P(*)^«(*) (C9) 
is also a solution given that: 

dso p{t) dp{t) 

— = , — — = -kso{t). (CIO) 

dt TO dt 

Now let us assume that the mass increases in time. We 
will use again the hydrodynamic equations (jC4p and 
(|C5|) . however contrary to the stationary problem dis- 
cussed above, the hydrodynamic approximation is a cru- 
cial assumption for having the center of mass mode in- 
dependent of the interaction strength. Strictly speaking, 
the Galilean invariance ljC9|) is valid only if mass, cur- 
vature and interaction remain constant in time. In the 
hydrodynamic regime the shape of the condensate does 
not depend on the mass, therefore we may seek a solution 
in a form similar to (|C7p : 



Pit) 



p-u^-m^ dp 



-fit)- 



(Cll) 



U dt 

Substituting these formulas into (|C4|I and IjCSfl yields 



di^ Mtj^ 



(C12) 
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which is to be suppHed with the initial conditions: p(0) = 
0, p{0) — —Xq k. In the adiabatic hmit this equation has 
an approximate solution: 

sin / uj{T)dT, (C13) 
Jo 

f{t) w xok ('^] ' cos /* u{T)dT. (C14) 



with ti'(T) = \J ■ The first of these equations shows 

that as the mass increases (or equivalently the tunneling 
constant decreases) the momentum also increases. 
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